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13.  This  report  presents  kinetics  models  for  the  chemically  facilitated 
gasification  of  the  ubiguitous  boron  oxide  coating  that  inhibits  particulate 
boron  ignition  and  the  second  stage  high  temperature  boron  surface  burning 
that  occurs  once  the  oxide  coating  has  been  removed.  These  models  include  a 
homogeneous  gas  phase  oxidation  mechanism  consisting  of  19  chemical  species 
and  58  forward  and  reverse  elementary  reactions,  multicomponent  gas  phase 
diffusion  and  heterogeneous  reactions  between  gas  phase  oxidation  products  and 
the  boron  oxide/boron  surface.  Kinetic  mechanisms  are  formulated  to  describe 
the  heterogeneous  surface  reactions  for  both  8203(1)  and  B(s)  surfaces  in 
terms  of  simple  adsorption  and  desorption  reaction  steps.  Thermochemical 
parameters  for  these  reaction  steps  are  estimated  from  the  current  heats  of 
formation  and  gas  phase  bond  energies  for  the  major  reactants.  Reaction  rate 
parameters  are  estimated  from  available  experimental  data  or  using  elementary 
transition  state  theory  arguments . 

The  kinetic  processes  are  treated  using  a  generalized  kinetics  computer  code, 
with  embedded  gradient  sensitivity  analysis,  for  chemically  reacting  flows. 
This  code  was  developed  to  provide  a  quantitative  model  for  the 
nonsteady-state  combustion  of  a  one-dimensional  (particle  radius)  spherical 
particle.  This  code  uses  a  non-linear,  implicit  partial  differential  equation 
solver  with  fully  automatic  time  step/mode  control  that  is  compatible  with 
formal  gradient  sensitivity  analysis  techniques. 
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(2)  Thermochemical  data  (entropies,  heats  of  formation  and  bond 
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(3)  A  one  dimensional,  nonsteady-state  model  for  the  combustion  of  a 
spherical  particle  was  formulated  and  the  associated  computer 
codes  developed.  This  model  can  currently  treat  particle 
shrinkage,  density  changes,  multi-component  gas  phase  diffusion  to 
the  particle  surface,  Stefan  flux  radiative  heat  transfer,  and 
homogeneous  (gas/gas)  and  heterogeneous  (gas/solid)  chemical 
reactions. 

(4)  Gradient  sensitivity  analysis  techniques  were  extended  to 
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and  reaction  enthalpies  and  rate  parameters  were  generated  from 
available  experimental  data  or  estimated  from  elementary 
transition  state  arguments. 

(6)  Model  results  were  obtained  for  the  gasification  of  a  liquid  boron 
oxide  particle  and  for  the  surface  oxidation  of  a  solid  boron 
particle. 


ABSTRACT 


This  document  is  the  final  report  to  the  Air  Force  Office  of  Scientific 
Research  on  work  performed  under  contract  F49620-88-C-0048.  It  covers  the 
period  between  January  1,  1988  and  December  31,  1990.  The  project's  primary 
objectives  were  to  model  hydrocarbon  assisted  boron  combustion  for 
air-breathing  systems  and  apply  gradient  sensitivity  analysis  to  determine  key 
reaction  pathways  and  critical  uncertainties  in  thermochemical  and  kinetic 
parameters. 

This  report  presents  kinetic  models  for  the  chemically  facilitated 
gasification  of  the  ubiguitous  boron  oxide  coating  that  inhibits  particulate 
boron  ignition  and  the  second  stage  high  temperature  boron  surface  burning 
that  occurs  once  the  oxide  coating  has  been  removed.  These  models  include  a 
homogeneous  gas  phase  oxidation  mechanism  consisting  of  19  chemical  species 
and  58  forward  and  reverse  elementary  reactions,  multicomponent  gas  phase 
diffusion  and  heterogeneous  reactions  between  gas  phase  oxidation  products  and 
the  boron  oxide/boron  surface.  Kinetic  mechanisms  are  formulated  to  describe 
the  heterogeneous  surface  reactions  for  both  6203(1)  and  B(s)  surfaces  in 
terms  of  simple  adsorption  and  desorption  reaction  steps.  Thermochemical 
parameters  for  these  reaction  steps  are  estimated  from  the  current  heats  of 
formation  and  gas  phase  bond  energies  for  the  major  reactants.  Reaction  rate 
parameters  are  estimated  from  available  experimental  data  or  using  elementary 
transition  state  theory  arguments. 

The  kinetic  processes  are  treated  using  a  generalized  kinetics  computer 
code,  with  embedded  gradient  sensitivity  analysis,  for  chemically  reacting 
flows.  This  code  was  developed  to  provide  a  quantitative  model  for  the 
nonsteady-state  combustion  of  a  one-dimensional  (particle  radius)  spherical 
particle.  This  code  uses  a  non-linear,  implicit  partial  differential  equation 
solver  with  fully  automatic  time  step/node  control  that  is  compatible  with 
formal  gradient  sensitivity  analysis  techniques. 
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1.0  INTRODUCTION 


1. 1  Statement  of  Problem  and  Research  Obiectives 

For  several  decades  elemental  boron  has  long  been  proposed  as  an  advanced 
fuel  for  combustion  and  propulsion  systems  because  of  its  high  energy 
density. Based  on  the  heats  of  combustion,  the  gravimetric  heating 
value  of  elemental  boron  is  approximately  35%  greater  than  JP4,  a  widely  used 
conventional  liquid  hydrocarbon  fuel,  and  approximately  85%  greater  than 
metallic  aluminum,  another  advanced  fuel  candidate  currently  used  in  solid 
fueled  rockets.  On  a  volumetric  basis,  boron  is  even  more  appealing  with  a 
heating  value  approximately  300%  greater  than  JP4  and  60%  greater  than 
aluminum.  These  potential  advantages  have  sparked  a  great  deal  of  interest  in 
the  enhancement  of  conventional  liquid  hydrocarbon  combustion  via  the  addition 
of  particulate  boron  in  the  form  of  slurry  suspensions.  For  air-breathing 
systems,  however,  the  theoretical  potential  has  not  been  fully  realized  in 
either  research  or  practical  combustors . 

Two  major  unresolved  combustion  problems  are  the  poor  combustion 
efficiency  and  inadequate  ignition  characteristics  of  boron  particles. ^“3  x0 
achieve  the  full  theoretical  performance  of  boron  combustion,  boron  particles 
must  ignite  and  completely  burn  (i.e.  reach  final  products  of  gaseous  H2O, 
gaseous  CO2  and  liquid  boria  B2O3)  within  a  very  limited  residence  time 
(2  -  10  ms).  The  poor  combustion  efficiency  results  from  slow  kinetically 
controlled  condensation  processes  to  form  the  liquid  phase  metal  oxide  which 
do  not  reach  the  desired  equilibria  during  the  short  residence  times  available 
in  most  combustion  devices.  The  ignition  problem  results  from  the  formation 
of  an  oxide  coating  on  unburned  boron  fuel  particles  which  inhibits  boron 
surface  oxidation  and  vaporization  during  the  initial  heat-up  phase  of  the 
particle. 
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Finding  adequate  solutions  to  these  practical  problems  has  been  hindered 
by  the  absence  of  a  fundamental  understanding  of  the  combustion  processes. 

The  current  understanding  of  boron  combustion  consists  of,  at  best,  an  overall 
description  of  the  most  probable  physical  and  chemical  processes;  i.e. 
knowledge  of  elementary  processes  and  parameters  is  significantly  lacking.  It 
is  apparent  from  previous  combustion  research,  that  if  solutions  to  the 
current  problems  preventing  practical  application  of  boron  combustion  are  to 
be  found,  additional  research  will  be  required,  particularly  at  the  elementary 
level. 

Recognition  of  the  research  needed  to  help  resolve  the  practical 
difficulties  associated  with  boron  combustion  has  motivated  several  on-going 
research  efforts  designed  to  obtain  a  better  understanding  of  the  combustion 
process  at  the  molecular  level.  These  include  experimental  studies  of  small 
boron  cluster  structure,  stability  and  reactivity, 1-4-1- 7  detailed  chemical 
kinetic  combustion  modeling, 1-8-1-15  flow  reactor  rate  measurements  for  bulk 
boron  oxidation  and  gasification  by  0(g),  C02(g),  8203(g)  and  H20(g) 1-16-1-21 
and  boron  surface  structure  and  reactivity  studies . 1“22, 1-23 

As  part  of  this  ongoing  research  effort,  the  objective  of  the  present 
program  was  to  construct  a  detailed  model  for  particulate  boron  assisted 
hydrocarbon  combustion.  The  goal  of  the  program  was  to  obtain  a  better 
fundamental  understanding  of  the  elementary  processes  and  the  role  these 
processes  play  during  the  entire  combustion  process.  Since  the  long  range 
success  of  this  program  would  be  dependent  on  the  ability  to  obtain  the 
necessary  data  for  input  to  such  a  model,  the  program  was  designed  to  be 
closely  coordinated  with  concurrent  experimental  work  by  other  groups.  In 
addition,  the  program  relied  on  gradient  sensitivity  analysis 

techniques ^”^4, 1-25  to  the  development  and  validation  of  the  model  and  to 

efficiently  interface  modeling  and  experimental  programs. 
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1.2  Program  Scope 


Boron-assisted  hydrocarbon  combustion  is  complex  process  which  includes  a 
heterogeneous  (solid/gas  and  solid/ liquid)  ignition,  volatization  and 
oxidation  phase,  a  homogeneous  (gas/gas)  oxidation  phase,  and  a  heterogeneous 
(gas/liquid)  condensation  phase  as  combustion  products  cool.l"^  jn  earlier 
work,  the  homogeneous  boron  combustion  phase  was  modeled  and  gradient 
sensitivity  analysis  was  used  to  determine  key  reaction  pathways  and  critical 
uncertainties  in  thermodynamic  and  kinetic  parameters.  1“12  The  primary 
goal  of  the  current  research  was  to  develop  a  combustion  model  for  the  initial 
heterogeneous  ignition,  volatization  and  oxidation  phase  of  boron  combustion 
and  to  couple  that  model  to  the  homogeneous  oxidation  model  developed  earlier. 
As  in  the  previous  work,  gradient  sensitivity  analysis  techniques 1“24, 1-25 
were  used  to  identify  key  reaction  mechanisms  and  critical  uncertainties  in 
molecular  thermodynamic  and  kinetic  parameters. 

1.3  Program  Overview 

Experimental  observations  indicate  that  the  combustion  of  boron  particles 
occurs  in  two  stages.  Macek,l“26  for  example,  observed  that  after  heating  to 
between  1800  K  and  2000  K,  the  boron  particle  becomes  luminous  and  glows  for  a 
short  period  of  time,  but  may  fade  out  as  the  boron  particle  is  weighted  with 
a  molten  boron  oxide  layer.  If  the  ambient  temperature  is  sufficiently  high, 
then  this  initial  ignition  state  is  followed  by  second  stage  of  burning  which 
is  much  brighter  and  longer  than  the  first  stage.  The  second  stage 
corresponds  to  the  full-fledged  combustion  of  a  "clean"  boron  particle  which 
occurs  once  the  oxide  coating  has  been  removed  by  vaporization  and 
heterogeneous  surface  reactions. 

Figure  1.1  depicts  schematically  the  two  stages  characterizing  the 
combustion  of  a  single  boron  particle  in  chemically  reactive  environments .  As 
indicated  in  this  figure,  a  detailed  simulation  of  the  combustion  process 
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requires  homogeneous  gas  phase  oxidation  reactions  and  multicomponent  gas 
phase  diffusion  to  the  particle  surface,  heterogeneous  surface  reactions 
between  gas  phase  species  and  the  oxide  (stage  1)  or  solid  boron  (stage  2) 
surface,  liquid  oxide  phase  diffusion  and  chemistry,  heterogeneous  reactions 
at  the  boron  oxide-solid  boron  interface,  Stefan  flux  from  the  particle 
surface  and  radiative  heat  transfer. 

Modeling  the  full  combustion  process  is  clearly  an  inherently  formidable 
task.  Consequently,  the  approach  of  the  current  program  was  to  develop  a 
complete  combustion  model  in  successive  stages  which  would  progressively  build 
on  each  other  by  incorporating  additional  important  processes.  Each  distinct 
stage  was  designed  to  focus  on  key  issues  concerning  the  combustion  process 
and  was  iteratively  upgraded  to  reflect  comparisons  with  phenomenological 
observations  and  new  experimental  or  theoretical  measurements  of 
thermochemical  and  kinetic  parameters. 

Accordingly,  a  simplified  description  of  boron  combustion  was  formulated 
which  focused  on  the  gasification  of  molten  boron  oxide  layer  and  the  second 
stage,  high  temperature  surface  burning  of  a  solid  boron  particle.  A 
schematic  presentation  of  these  sub-models  is  shown  in  Figure  1.2.  The  liquid 
oxide  droplet  model  is  a  first-order  representative  of  that  phase  of  boron 
combustion  when  the  boron  particle  is  coated  with  an  appreciable  oxide  layer. 
The  oxide  layer  is  removed  at  the  surface  through  vaporization.  However,  at 
high  temperatures  the  gasification  process  can  involve  chemical  transforma¬ 
tions.  Several  earlier  investigations ^-"3, 1-27  have  indicated  that  water,  for 
example,  can  promote  ignition.  Similarly,  because  of  its  low  volatility, 
gasification  of  particulate  boron  involves  chemical  transformation  at  the 
surface  into  more  volatile  suboxides  and  suboxyhydrides .  Thus,  one  key 
objective  of  this  submodel  was  to  reliably  predict  the  heterogeneous  kinetic 
mechanisms  by  which  surrounding  gases  significantly  after  the  rates  of  boron 
and  boron  oxide  gasification  evaporation. 
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During  the  first  year  of  this  program,  the  two  submodels  depicted  in 
Figure  1.2  were  treated  using  the  continuously  stirred  reactor  model 
(CSRM).1-13  This  is  a  steady-state  model  in  which  spatial  degrees-of-freedom 
are  eliminated  by  assuming  a  thin  control  volume  perpendicular  to  the  flow 
direction  which  is  well  mixed  so  that  species  profiles  are  uniform.  The  goal 
of  that  work  was  to  upgrade  the  homogeneous  gas  phase  oxidation  model  by 
incorporating  heat  and  mass  transport.  Mass  transport  was  treated  by  both 
(1)  parameterizing  the  mass  fluxes  by  the  residence  time  and  (2)  using  molar 
production  rates  for  heterogeneous  surface  reactions.  The  results  of  that 
analysis  were  presented  in  the  first  annual  report  and  will  not  be  reproduced 
here.1”13  The  principle  result,  however,  was  that  the  kinetic  calculations, 
coupled  with  sensitivity  analysis,  indicated  that  the  reactions  identified 
previously1”11* 1-12  as  being  important  in  the  homogeneous  oxidation  model  for 
a  closed  system  remain  the  dominant  mechanisms  when  heat  and  mass  fluxes  are 
inc luded. 

During  the  second  year,  the  CSRM  was  superseded  with  a  detailed 
one-dimensional,  time  dependent  spherical  particle  combustion  model  and  the 
associated  computer  codes  which  include  the  capability  to  compute  sensitivity 
coefficients  were  completed.1-1^  This  model  represents  state-of-the-art 
modeling  technology  and  marks  a  significant  achievement  in  the  proposed  goals 
of  this  project.  The  spherical  particle  model  can  describe  particle 
shrinkage,  density  changes,  Soret  effects  and  radiative  heat  transfer.  In 
addition,  the  ability  to  treat  surface  speciation  and  elementary  absorption 
and  desorption  reactions  represents  an  unsurpassed  level  of  detail  in  modeling 
such  complex  systems.  As  in  earlier  work,  gradient  sensitivity  analysis 
techniques  were  incorporated  into  the  computer  codes  to  facilitate  the 
identification  of  key  mechanisms  and  system  parameters,  as  well  as  to  provide 
an  internal  mechanism  to  aid  in  evaluating  the  validity  and  self-consistency 
of  the  model  input  parameters.  This  extension  of  gradient  sensitivity 
analysis  techniques  to  heterogeneous  processes  such  as  particle  ignition  was 
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not  only  critical  to  the  successful  treatment  of  boron  combustion,  but  will  be 
a  major  contribution  to  models  of  other  inherently  heterogeneous  chemical 
systems . 

During  the  third  year,  the  most  potentially  important  heterogeneous 
surface  reactions  were  identified,  first-order  formulations  of  the  reaction 
mechanisms  were  developed  and  thermochemical  and  reaction  rate  parameters  were 
estimated.  These  surface  reaction  mechanisms  include  adsorption  and 
desorption  processes  and,  therefore,  differ  considerably  from  earlier  modeling 
work  in  their  level  of  microscopic  detail.  This  was  motivated  by  the  desire 
to  ultimately  provide  a  more  realistic  model  of  boron  combustion  which,  for 
example,  would  provide  a  description  of  surface  speciation  and  a  more  accurate 
dependence  of  reaction  rates  on  species  concentration.  Thus,  while  the 
proposed  mechanisms  are  presently  very  speculative  because  of  the  lack  of 
detailed  experimental  or  theoretical  data  for  elementary  processes,  they 
provide  the  needed  first  step  in  developing  an  accurate  combustion  model.  In 
addition,  the  detailed  level  of  description  complements  concurrent 
experimental  studies  on  small  boron  c lusters 1“4-1~7  by  Anderson's  group  at 
Stony  Brook  and  boron  surface  structure  and  reactively1-^,  1-23  by  Trenary's 
group  at  Chicago.  Thus,  the  reaction  mechanisms  presented  in  this  report 
should  be  viewed  as  an  initial  formulation  which  can  be  modified  and 
interatively  upgraded  as  new  experimental  data  becomes  available. 

l.A  RspsdL-SLvt.Une 

Section  2.0  presents  the  mathematical  formulation  of  the  combustion 
process.  This  includes  the  governing  equations  for  mass  and  energy 
conservation,  as  well  as  a  description  of  the  numerical  techniques  used  to 
solve  these  equations.  The  formalism  and  computational  methods  for  applying 
gradient  sensitivity  analysis  techniques  to  heterogeneous  processes  with 
spatial  and  temporal  dependence  are  also  presented  in  Section  2.0. 
Heterogeneous  surface  reactions,  reaction  mechanisms  and  estimates  of 
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thermochemical  and  parameters  are  presented  in  Section  3.0.  Section  4.0 
presents  model  results  for  boron  oxide  gasification  and  particulate  boron 
surface  burning  in  hydrocarbon  combustion  environments. 
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2.0  ONE  DIMENSIONAL  SPHERICAL  PARTICLE  COMBUSTION  MODEL 


This  section  describes  the  mathematical  formulation  of  a  one-dimensional, 
nonsteady- st ate  model  for  the  combustion  of  a  spherical  particle. 

Sections  2.1  -  2.3  present  mass,  species  and  energy  conservation  equations  and 
Section  2.4  specifies  the  appropriate  boundary  conditions.  In  Section  2.5, 
the  numerical  methods  used  to  solve  these  strongly  coupled  partial 
differential  equations  (PDE)  are  discussed  and  qualitatively  compared  with 
other  frequently  used  numerical  techniques.  The  application  of  sensitivity 
analysis  techniques  to  problems  in  which  system  variables  and  parameters  are 
both  temporally  and  spatially  dependent  is  discussed  in  Section  2.6  and 
includes  a  description  of  the  computational  method  used  to  calculate 
sensitivity  coefficients. 

2.1  Governing  Equations  for  the  Gas  Phase 

The  governing  equations  for  overall  mass,  species,  and  energy 
conservation  in  the  gas  phase  are 
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where  a'  represents  a  geometrical  factor  (a  =  0  for  cartesian  coordinates, 
a  -  1  for  cylindrical  coordinates,  and  a  »  2  for  spherical  coordinates),  pg 
denotes  the  density  of  the  gas  mixture,  Cpjg  the  specific  heat  capacity  of  the 
gas  mixture  at  constant  pressure,  Xg  the  thermal  conductivity  of  the  gas 
mixture,  Ygji  the  mass  fraction  of  the  i**1  species  in  the  gas  phase,  Tg  the 
temperature  of  the  gas  mixture,  Cp^i  the  specific  heat  capacity  of  the  ith 
gas  phase  species,  Hgi  the  specific  enthalpy  of  the  ith  species  with  respect 
to  the  mixture  of  remaining  species,  Vri  the  diffusion  velocity  of  the  ith 
species,  vr  the  fluid  velocity  in  the  radial  direction,  r  the  radial 
coordinate,  and  Rgsi  the  mass  production  rate  of  the  i**1  species  by  gas  phase 
chemical  reaction.  The  ideal  gas  equation  of  state  relates  Y„  j,  p  and  T„  by 
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where  P  is  the  thermodynamic  pressure,  R°  the  universal  gas  constant,  and  W £ 
the  molecular  weight  of  the  ith  species. 
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The  transport  model  recommended  by  Coffee  and  Heimerl^"!  is  used  to 
describe  multi-component  diffusion  processes.  According  to  Coffee  and 
Heimerl^"!  the  diffusive  flux  is  described  by 
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where  dr>j[  is  the  ordinary  diffusion  velocity,  is  the  thermal  diffusion 

velocity,  and  dc  is  the  constant  diffusion  velocity.  The  ordinary  diffusion 
coefficient  can  be  conveniently  calculated  by  using  an  effective  binary 
diffusivity  as  follows. 
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where  is  the  mole  fraction  of  the  i**1  species  in  the  gas  phase  and  is 

dependent  on  radial  position.  In  order  to  generalize  the  binary  diffusion 
formulas  and  mass  transfer  correlations  by  simply  using  the  Dj,  it  is  assumed 
(Curtiss  -  Hirshfelder  approximation)  that  all  but  the  i^  species  move  with 
the  same  velocity.  Then  the  Dj  can  be  related  to  the  binary  diffusion 
coefficients  Djf  through  the  following  relation. 
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The  thermal  diffusion  velocity,  dt>i,  is  obtained  by  employing  the  thermal 
diffusion  ratio  (tp^). 
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The  constant  diffusion  velocity,  dc,  is  introduced  to  assure  the  zero  net 
diffusive  flux. 

j=n 
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The  transport  coefficients  are  calculated  by  using  the  Transport  Package 
developed  by  Kee  et  al. ,2~2  which  evaluates  gas  phase  conductivities  and 
diffusion  coefficients  from  the  Lennard-Jones  potential,  the  well  depth,  the 
collision  diameter,  the  dipole  moment,  the  polarizability,  and  the  rotational 
relaxation  collision  number. 

The  inner  boundary  located  at  the  droplet  surface  moves  as  the  droplet 
surface  regresses.  Thus,  a  density  weighted  coordinate  transformation  is 
introduced  to  immobilize  the  regressing  interface  location.  Specifically, 
defining  <f>  by 
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The  gas  phase  density  (pg)  is  usually  much  smaller  than  the  condensed  phase 
density  (pc),  i.e.  pg  «  pCj  except  in  the  critical  conditions,  ~_.:h  that  the 
mp  becomes  identical  to  the  mass  flow  rate.  The  following  relationships 
derived  from  the  chain  rule  are  used  to  perform  the  coordinate  transformation: 
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The  continuity  equation.  Equation  (2-1),  is  automatically  satisfied  in  the  new 
coordinate  system.  Using  Equations  (2-14)  and  (2-15),  the  species  and  energy 
conservation  equations.  Equations  (2-2)  and  (2-3),  become 
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2 . 2  Governing  Equations  for  the  Condensed  Phase 

The  size  of  the  condensed  phase  (oxide  droplet  or  boron  particle)  is 
continuously  decreasing  with  time  as  combustion  or  vaporization  proceeds.  The 
overall  mass  conservation  condition  determines  the  rate  of  change  of  the 
particle  diameter  by 
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where  pc  is  the  particle  density.  The  energy  conservation  equation  becomes 
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where  the  subscript  c  denotes  the  condensed  phase,  rp  the  particle  radius, 
Cp>c  the  specific  heat  capacity  and  Xc  the  thermal  conductivity. 


3T 

p  c  r-* 

c  p,c  3t 


a-1  a-1 

"p  n 


£  T 

3_  a-1  lie,  _  In 

3ri  c11  3ii  ^  pc  p,c  3t 


3T 

_ £ 

3q 


(2-20) 


2-6 


2.3  Governing  Equations  for  the  Surface  Species 


In  the  present  surface  chemistry  model,  chemical  reactions  among  surface 
species  or  between  surface  species  and  gaseous  species  are  not  considered. 
Therefore,  the  governing  equations  for  the  surface  species  can  be  written  as 
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where  Vg^j^  denotes  the  stoichiometric  coefficient  for  species  i  appearing  as 
a  adsorption  product  of  species  k,  ka>k  the  adsorption  rate  constant  for 
gaseous  species  k,  k^i  the  desorption  rate  constant  for  surface  species  i, 
CSji  the  concentration  of  surface  species  i. 

2.4  Initial.  Boundary,  and  Interface  Conditions 

In  order  to  solve  the  governing  equations  described  above,  the  attending 
initial,  boundary,  and  interface  conditions  should  be  specified.  Initial 
conditions  specifies  the  species  concentration  and  ambient  and  particle 
temperatures  at  t  =  0,  i.e. 
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(2-25) 


Tc(0,r)  =  Tc>o(r)  0  <  r  <  rp 

The  boundary  conditions  at  the  left  hand  side,  which  is  positioned  at  the 
center  of  the  condensed  phase,  are  based  on  no  flux  conditions  for  both 
chemical  species  and  temperature  of  the  condensed  phase.  Thus,  at  r  =  0 

3T 

7-^=0  (2-26) 

3r 

The  boundary  conditions  at  the  right  hand  side  are  assumed  to  be  same  as  the 
ambient  conditions  so  that  at  r  =  « 

Y  =  Y  (2-27) 
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The  temperature  profile  is  continuous  across  the  interface, 

at  r  =  r  ,  T  =  T  .  (2-29) 

P  g  s 

The  interface  conditions  at  the  surface  of  the  condensed  phase  are  determined 
by  considering  the  mass  balance  across  the  boundary;  the  chemical  species  flux 
at  the  surface  is  same  as  the  amount  destroyed  or  produced  by  the  surface 
chemical  reactions 


2-8 


(2-30) 
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Summing  the  above  equations  over  chemical  species  and  using  the  fact  that 
summation  of  diffuse  fluxes  should  be  equal  to  zero,  the  mass  flux  (pg  vr)  at 
the  surface  can  be  written. 
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Similarly,  the  energy  balance  across  the  interface  yields 
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2.5  Numerical  Methods 

In  Sections  2.1  thru  2.4,  the  vaporization  and  combustion  of  isolated 
droplets  are  treated  as  an  initial  boundary  value  problem  in  which  governing 
parabolic  PDEs  are  strongly  coupled  due  to  strong  chemical  interactions. 
General  solution  procedures  for  initial  boundary  value  problems  are,  first,  to 
discretize  the  spatial  part  of  the  differential  equations  and,  then,  to 
integrate  the  resulting  system  of  ordinary  differential  equations  over  time. 

2.5.1  Temporal  Integration 

Widely  used  numerical  techniques  for  temporal  integration  in  chemically 
reactive  flow  problems  include  time  splitting  techniques, 2”3  semi-implicit 


2-9 


methods, ,2-5  Newton  linearization  techniques,  and  Gear  type  non-linear 
iterations . 2“6 

A  high  computational  efficiency  may  be  achieved  by  using  the  time 
splitting  technique  in  solving  chemically  reactive  flow  problems. 2-7-2- 10  The 
time  splitting  technique  breaks  complex  interactions  among  processes  into 
individual  processes  such  that  each  individual  process,  especially  chemical 
kinetics,  can  be  solved  separately.  Therefore,  the  number  of  simultaneous 
equations  to  be  solved  can  be  dramatically  reduced.  For  problems  with 
disparate  time  scales  among  sub-processes,  an  efficient  temporal  integration 
is  achieved  by  applying  the  appropriate  time  step  according  to  the 
characteristic  time  scale  of  each  sub-process.  Thus,  for  chemically  reactive 
flow  problems  in  which  chemical  processes  are  much  faster  than  transport 
processes,  a  small  time  step  is  applied  to  resolve  chemistry  equations, 
whereas  transport  equations  can  be  integrated  with  a  large  time  step.  In 
addition,  a  high  degree  of  parallelism  can  be  achieved  because  chemistry 
equations  at  each  grid  point  are  decoupled  and  solved  in  parallel.  This  makes 
the  time  splitting  techniques  extremely  favorable  to  the  parallel  super 
computer  architectures  such  as  the  Connection  Machine  and  the  Massively 
Parallel  processes. 2-11  However,  numerical  accuracy  may  be  considerably 
degraded  when  solutions  from  individually  splitted  differential  equations  are 
combined  to  yield  real  solutions.  Furthermore,  as  the  system  approaches 
extinction,  the  time  scale  of  transport  processes  becomes  comparable  to  or 
even  smaller  than  that  of  chemical  processes,  and  the  interaction  between 
transport  and  chemical  processes  become  important.  Thus,  the  time  splitting 
technique  which  decouples  these  processes  may  not  resolve  the  transition  from 
quasi-steady  burning  to  extinction  accurately. 

A  Newton  linearization  technique2-12-2-13  or  semi-implicit  method^-^  may 
be  applied  to  decouple  highly  non-linear  equations  appearing  in  chemically 
reactive  flow  problems.  Numerical  comparisons  between  these  two  techniques  in 
solving  transport/reaction  equations  were  made  by  Yuk.2"14  These  techniques 
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may  not  resolve  numerical  stiffness  problems  arising  from  chemical  species 
with  high  reaction  fluxes  such  as  radicals.  Because  the  chemical  species  with 
high  reaction  fluxes  are  very  likely  to  be  in  the  steady  state,  these 
numerical  problems  may  be  avoided  by  carefully  applying  the  pseudo-steady 
state  approximation  to  remove  numerically  stiff  differential  equations  from 
the  ODE  systems.  However,  if  used  inappropriately,  the  steady  state 
approximation  may  be  a  source  of  a  significant  error. 

In  the  present  study,  Gear's  implicit  method^-**  is  employed  to  integrate 
systems  of  stiff  ODEs  over  time.  Gear's  method  automatically  controls  the 
time  step  and  the  order  of  integration  to  maintain  requisite  levels  of 
stability  and  accuracy  in  the  time  domain.  Even  though  Gear's  method  yields 
accurate  numerical  solutions  compared  to  the  other  techniques  described  above, 
it  requires  inversion  of  a  large  matrix,  which  is  a  computationally  intensive 
operation.  It  is  hoped  that  the  numerical  code  developed  here  using  Gear's 
method  may  be  used  to  critically  evaluate  the  code  with  a  more  computationally 
efficient  temporal  integrator  such  as  the  time  splitting  technique,  the 
semi-implicit  method,  and  the  Newton  linearization  in  solving  chemically 
reactive  flow  problems. 

2.5.2  Spatial  Discretization 

The  governing  PDEs  given  in  Sections  2.1  thru  2.3  are  spatially 
discretized  by  a  Finite  Element  method.  The  particular  Finite  Element  model 
used  in  this  study  employs  a  Galerkin  formulation  and  linear  basis  functions. 
Unlike  conventional  Finite  Difference  methods,  the  flux  boundary  conditions 
are  exactly  incorporated  into  the  Finite  Element  formulation  without 
introducing  any  fictitious  model  nodal  points  or  degrading  the  numerical 
accuracy.  This  allows  accurate  calculation  of  flux  boundary  conditions,  which 
is  very  important  especially  in  the  present  case  because  it  leads  to  an 
accurate  prediction  of  the  particle  gasification  rate. 
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In  chemically  reactive  flow  problems  including  the  vaporization  and 
combustion  of  isolated  particles,  one  or  more  steep  fronts  may 
appear/disappear  or  move  continuously  in  the  solution  domain.  The  spatial 
stiffness  may  be  resolved  by  increasing  the  order  of  basis  functions  in 
Orthogonal  Collocation  or  Galerkin  Finite  Element  formulations  in  the  vicinity 
of  the  steep  fronts.  Even  though  further  studies  are  required  to  evaluate 
this  technique,  our  study  shows  that  when  Orthogonal  Collocation  Finite 
Element  method  is  implemented  with  adaptively  varying  the  degree  of  polynomial 
orders,  numerical  performance  is  very  good  in  terms  of  reducing  Gibbs 
overshooting/undershooting,  numerical  diffusion,  and  temporal  stiffness  due  to 
close  adjacent  grid  points.  This  technique  is  expected  to  perform  as  well 
with  the  Galerkin  Finite  Element  method.  However,  the  matrix  structure 
becomes  irregular,  which  causes  difficulties  in  computer  program  coding  and 
data  allocation,  and  considerably  increases  the  computational  time. 

The  alternate  approach,  which  is  chosen  in  the  present  study,  is  to  apply 
dense  gridding  in  the  vicinity  of  steep  fronts.  The  node  control  is  performed 
by  various  adaptive  gridding  techniques.  Philosophically,  adaptive  gridding 
techniques  widely  used  in  chemically  reactive  flow  problems  are  of  two  types. 
In  the  first  method,  adaptive  gridding  is  based  on  a  single  fluid  property  or 
a  single  dependent  variable  behavior.  This  is  implemented  by 

1.  Moving  nodes  at  mean  fluid  velocities  or  at  some  other 
characteristic  velocity  in  a  fluid,  which  is  known  as  Lagrangian 
regridding  procedures,  or 

2.  Performing  coordinate  transformation  to  equally  distribute  gradients 
of  important  variables  such  as  temperature. 

In  the  combustion  of  isolated  droplets,  Lagrangian  regridding  procedures  are 
inappropriate  because  neither  a  flame  velocity  nor  a  flame  position  is  well 
defined.  The  coordinate  transformation  based  on  the  equi-distribution  of 
gradients  of  temperature  was  tried  by  the  authors.  Severe  difficulties  were 
found  in  obtaining  convergence  in  the  Newton-Raphson  iteration  of  the 
non-linear  initial  value  problem  solver.  This  problem  may  be  avoided  by 
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applying  a  Newton  linearization  technique  to  perform  temporal  integration  as 
shown  by  Dwyer  and  Sanders^- but  such  a  temporal  integration  is  not 
acceptable  in  the  present  study  as  discussed  in  the  previous  section.  The 
poor  convergence  in  the  non-linear  initial  value  problem  solver  is  probably 
because  the  coordinate  transformation  may  increase  stiffness  or  because  the 
resulting  grid  system  which  is  based  only  on  the  temperature  profiles  may  not 
be  good  enough  to  resolve  all  the  steep  fronts  in  the  model  domain. 

The  second  method  considers  more  than  one  dependent  variable  to  perform 
adaptive  gridding.  This  is  implemented  by 

1.  Inserting  or  removing  grid  points  to  reduce  local  errors  or 
equi-distribute  mesh  functions,  or 

2.  Moving  nodes  smoothly  at  each  time  step  to  reduce  global  errors  or 
equi-distribute  mesh  functions. 

The  former  technique  (1)  is  based  on  two  step  procedures;  a  solution  step  and 
adaptive  gridding  step,  whereas  in  the  latter  technique  (2),  termed  moving 
node  technique,  the  node  positions  are  computed  together  with  the  solution 
values  at  each  node.  Thus,  node  control  is  performed  statically  in  the  former 
technique  (1),  and  dynamically  in  the  latter  technique  (2).  Adaptive  gridding 
based  on  inserting  or  removing  nodes  can  potentially  solve  partial 
differential  equations  to  a  prescribed  level  of  accuracy,  whereas  moving  node 
methods  can  follow  evolving  non-uniformities  and  very  effectively  reduce 
dispersive  errors  and  Gibbs  ov  rshooting/understooting.  These  adaptive 
gridding  techniques  are  different  from  adaptive  mesh  refinement 
techniques, 2-15-2-18  £n  which  appropriately  graded  meshes  are  constructed  from 
a  coarse  mesh  through  iterative  procedures.  In  general,  the  adaptive  gridding 
techniques  discussed  here  may  be  more  suitable  for  time  dependent  PDEs  than 
the  adaptive  mesh  refinement  techniques  because  they  modify  the  previous  time 
step  grid  system  slightly  at  each  time  step  according  to  evolution  of 
solutions. 
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Remarkable  progresses  have  been  made  toward  achieving  adaptive  gridding 
by  inserting  or  removing  grid  points  according  to  error  estimates .2-19  This 
adaptive  grid  technique  is  implemented  by  adding  or  removing  OOEs  to  the  model 
and  interpolating  or  projecting  dependent  variables  to  provide  initial 
conditions  for  the  added  ODEs.  Thus,  dependent  variable  interpolation  or 
projection  must  be  well  defined  and  performed  accurately  for  maximum 
effectiveness.  Whenever  nodes  are  added  or  deleted,  the  initial  value  problem 
solver  should  be  restarted,  which  is  very  costly  for  multi-step  methods. 

In  the  present  study,  moving  node  techniques  were  chosen  to  perform 
adaptive  gridding.  The  moving  node  techniques  move  a  fixed  number  of  nodes 
smoothly  at  each  time  step  to  adjust  nodes  following  evolution  of  solutions. 

No  dependent  variable  transformation  is  involved,  and  nodes  are  not  inserted 
or  removed.  Here,  the  moving  node  techniques  are  implemented  onto  the 
Galerkin  Finite  Element  method  to  formulate  a  Moving  node  Finite  Element 
Method  (hereafter,  M-FEM) .  For  convenience,  the  conventional  fixed  node 
Galerkin  Finite  Element  Method  is  termed  the  Fixed  node  Finite  Element  method 
(hereafter  F-FEM).  Unlike  the  F-FEM,  the  M-FEM  should  account  for  effects  of 
node  movements  on  dependent  variables.  The  successful  M-FEM  formulation 
should  be  able  to  control  nodes  in  a  non-stiff  manner,  and  prevent  any 
abnormal  behaviors  such  as  node  crossing.  A  detailed  mathematical  formulation 
is  shown  in  the  next  section. 

2.5.3  M-FEM  Formulation 

The  mathematical  derivations  presented  below  are  obtained  by  closely 
following  Miller  and  Miller's^^O  ^  Linch's2"21  work,  in  which  more  detailed 
discussions  with  rigorous  mathematical  justification  are  found. 
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In  general,  time  dependent  differential  equation  can  be  written  as 


u'  =  L  (u) 


t  >  0 


(2-33) 


where  L  (u)  is  a  non-linear  partial  differential  operator.  Because  nodes  are 
continuously  deforming  with  time,  the  approximate  solution,  v,  can  be  written 
in  terms  of  both  temporally  varying  node  positions,  s^(r),  and  nodal 
amplitudes,  a^j(t). 


v(t)  =  vCa1  (t),  ....a1  (t),  ...,an  (t),  ...an  (t),  s.U),  ...s  (t)) 
l  m  1  m  1  m 
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In  the  finite  element  r-ethod,  the  approximation  solution,  v,  is  given  by 

j=m 

vk(x,t)  *  J  ak  (t)  °j (x,t)  k  *  1,  ...,n  (2-35) 

j-1 

where  akj  is  a  nodal  amplitude  of  the  k**1  dependent  variable  at  the  j*-*1  node, 
m  and  n  are  the  number  of  node  points  and  dependent  variables,  respectively. 
cij(x,t)  is  the  corresponding  basis  function.  In  contrast  to  the  F-FEM,  the 
basis  function,  aj,  is  a  function  of  both  time  and  space  because  node  points 
move  continuously  with  time.  The  temporal  variation  of  v  is  derived  as 
follows  by  applying  the  chain  rule  differentiation  to  v  with  respect  to  time, 

.  k  j-m  9ak  as.  . 

if  -  l  “j  +  iT1  9  •  <2-36> 
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which  implies  that  the  approximate  solution  v  is  a  function  of  both  the 
amplitudes,  a^j ,  and  the  nodal  positions  sj.  The  basis  functions  txj  and 
are  defined  by 


a. 
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(2-37) 
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With  the  temporal  variation  of  v  given  in  Eq.  (2-36),  the  standard  method  of 
weighted  residual  method  is  applied  to  derive  the  M-FEM  formulation,  which  is 
given  by. 
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2.5.4  Nodal  Equations 

In  order  to  close  the  M-FEM  formulation,  it  is  necessary  to  describe  the 
temporal  variation  of  node  positions.  There  are  two  approaches  to  derive 
nodal  equations.  The  first  approach  is  to  equi-distribute  mesh  function 
values  over  all  zones. 2-22  The  second  approach  is  to  minimize  global 
errors . 2-20-2-23 
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In  the  mesh  function  approach,  mesh  functions  are  constructed  to  be  a 
measure  of  local  adequacies  of  mesh  deployment.  The  mesh  function  sought  is  a 
non-negative  and  dimensionless  function  with  the  following  properties, 

lim  8i  =  0  (2-40) 

h.— >0 

l 

where 


0*  «  1;  good  node  deployment 


0*  »  1;  poor  node  deployment 


The  node  moving  strategy  is  to  move  the  nodes  so  as  to  equi-distribute  the 
mesh  function  values  over  all  zones.  This  can  be  accomplished  by  solving  the 
following  equation. 
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Because  two  end  nodes  s^  and  sm  are  fixed,  i.e.,  ds^/dt  =  dsm/dt  =  0,  the 
above  differential  equation  can  be  easily  solved  to  describe  node  movements. 
The  mesh  function,  0-p1,  is  defined  as  linear  combination  of  appropriate 
sub-mesh  functions.  The  sub-mesh  functions  chosen  in  this  study  are 
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where  Qj  is  the  amount  of  heat  release  due  to  chemical  inactions  at  the  i1"*1 
node.  The  weights,  w^,  are  employed  to  control  the  influence  of  dependent 
variables  on  node  motion.  Each  is  chosen  to  normalize  the  mesh  function 
such  that 


i=m-l 

l  0*  =  1  (2-47) 
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The  first  and  second  mesh  functions,  which  were  originally  proposed  by  Madsen, 
are  designed  to  prevent  a  stiff  change  in  first  and  second  order  spatial 
derivatives.  The  third  mesh  function  is  implemented,  in  this  work,  to 
increase  grid  resolutions  in  the  vicinity  of  a  reaction  zone.  Then,  the  total 
mesh  function  becomes 

j=3 

0*  =  l  e.ej  (2-48) 

j=l 

where  ej's  are  user  specified  constants.  The  suitable  set  of  ej/s  is 
determined  only  by  trial-and-error  types  of  approaches.  The  values  of  ej 
chosen  in  this  study  are;  e^  =  1,  e£  =  0.005,  and  63  =  0.2.  Because  the  above 
mesh  functions  do  not  provide  resistances  against  node  crossing,  additional 
node  control  is  necessary  to  prevent  nodes  from  crossing  over  each  other  and 
maintain  a  good  node  topology.  Madsen  recommended  the  following  mesh  function 
to  control  over  excessive  adjacent  zone  distortion. 
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The  above  mesh  function  causes  the  nodes  to  be  equally  spaced.  However,  a 
non-uniform  node  spacing  may  be  more  suited  for  the  most  of  problems.  Thus, 
in  the  present  study,  the  additional  node  control  is  achieved  by  imposing  the 
following  constraint. 
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where  B  is  a  constant  larger  than  1.  The  B  is  chosen  to  be  5  in  this  study. 

It  is  important  to  note  here  that  the  values  chosen  for  e^'s  and  B  do  not 
affect  the  accuracy  of  model  solutions  significantly  but  do  influence  the 
numerical  efficiency.  In  addition  to  node  distortion,  it  should  be  noted  that 
a  extremely  high  node  density  in  a  small  region  may  increase  the  temporal 
stiffness  considerably.  Thus,  a  minimum  node  separation  distance  should  be 
enforced  to  maintain  the  numerical  stiffness  properties  of  the  ODE  system  at 
the  very  manageable  levels,  i.e. 


h.  <  h  . 

1  min 


(2-51) 


The  constraints  are  implemented  by  modifying  m^  in  small  and  distributing  the 
excess  over  zones.  The  actual  procedures  of  this  node  control  algorithm  are 
illustrated  by  Madsen. 2-22 

The  second  approach  is  proposed  by  Miller. 2-20,  2-24  This  technique 
implicitly  calculates  node  locations  to  minimize  global  errors.  Galinas  et 
al.2”25  derived  the  following  equations  by  minimizing  the  norm  of  the 
penalized  residual  function  with  respect  to  a^j  and  s^. 
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where  X^,  X^  are  penalty  functions  which  are  employed  to  prevent  the 
formulation  from  becoming  singular,  to  control  the  additional  numerical 
stiffness  arising  due  to  introduction  of  nodal  equations,  and  to  provide 
resistances  to  node  crossing.  Considerable  efforts  have  been  made  to  search 
for  appropriate  penalty  functions2”25-2-26  or  tQ  find  a  way  to  eliminate  the 
necessity  of  penalty  functions . 2-27 

The  performance  of  this  approach  with  various  penalty  functions  on  the 
problems  typically  arising  in  the  combustion  system  compares  very  well  with 
alternate  numerical  techniques2"28  solving  both  stable  and  unstable  freely 
propagating  premixed  laminar  flames  with  one  step  chemistry. 2-29  However,  as 
also  shown  by  Cho  et  al.2“29  for  the  simulation  of  droplet  combustion,  nodes 
are  often  locked  into  a  sharp  moving  front.  Thus  as  the  front  advances  across 
the  domain,  the  nodes  ahead  concentrate  and  the  nodes  behind  stretch 
drastically.  These  difficulties  may  be  resolved  in  a  certain  degree  by 
incorporating  Madsen's  mesh  function  approach  into  the  penalty  function. 
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(2-53) 


X2  = 
1 


i^m-1 


(As.  -  6)  +  C2  6XP  (  maX  (0‘0T)1) 
x  '  1=1 


X.w . 


l  i 


C2  i=m-l 

- -  +  C,  exp  (  max 

(ASi  -  6)  i=l 


(0-6^*)  (0-0*) 


(2-54) 


It  can  be  readily  shovm  that  these  penalty  functions  have  node  movements 
guided  by  the  mesh  function  approach,  when  the  node  deployment  violates 
seriously  the  ideal  mesh  system  defined  by  the  mesh  function. 

Numerical  experiments  have  been  performed  to  compare  the  above  two 
approaches.  It  was  found  that  the  second  approach,  which  is  based  on 
minimization  of  global  errors,  depends  too  much  on  the  selection  of  penalty 
functions  and  often  introduces  a  severe  numerical  stiffness.  On  the  other 
hand,  the  first  approach,  which  is  based  on  equi-distribution  of  mesh 
functions,  provides  only  a  sub-optimal  node  control,  but  no  additional 
numerical  stiffness  arises  from  the  node  control.  Because  the  problem  itself 
is  very  stiff  temporally,  any  additional  stiffness  due  to  a  node  control  may 
introduce  a  serious  difficulty  in  temporal  integration.  Thus,  the  first 
approach  is  employed  to  calculate  node  movements  in  the  present  study. 

2.5.5  Solution  Procedure 

The  combustion  and  vaporization  of  droplets  are  described  by  a  set  of 
differential-algebraic  equations,  Eq.  (2-10),  (2-16),  (2-17),  (2-18),  (2-20), 
(2-21),  and  (2-31).  The  gas  phase  mass  and  heat  balance  equations.  Equations 
(2-l<  j  and  (2-17)  are  descretized  by  the  M-FEM.  The  liquid  phase  equations, 
Eq.  (2-21)  are  descretized  by  the  F-FEM  with  uniformly  spaced  gridding.  Then, 
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the  resulting  ODE  systems  are  integrated  by  coupling  with  the  overall  mass 
balance  equation  of  the  liquid  droplet,  Eq.  (2-18),  and  the  concentrations  of 
the  surface  species. 

In  addition  to  chemical  species  concentration  and  temperature  at  each 
grid  point  and  the  droplet  radius,  mp  and  r  are  unknown.  Among  unknowns,  the 
mp  is  a  single  unknown,  but  appears  in  almost  every  equation.  In  order  to 
maintain  the  banded  matrix  structure,  the  mp  is  treated  as  unknown  at  each 
grid  point.  Then,  a  set  of  algebraic  equations  is  added  into  the  ODE  systems 
to  assure  that  these  variables  are  equal  everywhere.  A  similar  approach  was 
taken  by  Smooke  et  al.^“16  to  determine  the  flame  speed  in  freely  propagating 
adiabatic  flame.  The  resulting  ODE  system  with  the  tridiagonal  matrix 
structure  is  solved  by  the  various  implementations  of  Gear's  method. 

In  the  most  previous  approaches,  the  Jacobian  matrix  is  evaluated 
numerically  to  reduce  efforts  for  computer  program  coding.  In  the  current 
work,  the  analytical  form  of  Jacobian  matrix  is  used  to  improve  convergeance 
as  well  as  a  numerical  convergeance  for  sensitivity  coefficients. 

Furthermore,  it  was  found  that  the  usage  of  the  analytical  Jacobian  matrix 
reduces  computational  time  for  Jacobian  matrix  evaluation  significantly.  As  a 
result,  the  total  computation  time  required  for  the  solutions  was  reduced  by 
over  10  times. 

Even  though  all  the  results  presented  in  this  paper  are  based  on  the 
above  procedure,  it  is  worthwhile  to  note  that  the  energy  conservation 
equation  for  the  solid  phase  may  be  decoupled  from  those  for  the  gas  phase 
without  loss  of  significant  accuracy,  because  the  liquid  phase  processes  are 
much  slower  than  the  gas  phase  processes. 
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2.6  Sensitivity  Analysis 


Sensitivity  analysis  provides  a  formal  method  for  extracting  the 
relationship  between  input  and  output  parameters.  Sensitivity  analysis  can 
be  used  to  evaluate  the  sensitivity  of  output  parameters  to  model  input 
parameters  as  well  as  the  uncertainty  in  their  values.  Hence,  the  sensitivity 
analyses  would  be  powerful  interpretive  tools  in  investigating  chemically 
reactive  flow  which  is  complicated  by  non-linear  interactions  among  various 
processes  and  includes  a  large  number  of  input  parameters  containing 
significant  uncertainties.  However,  the  application  of  the  sensitivity 
analyses  especially  to  the  space-time  problem  may  become  inhibitely  expensive 
greatly  exceeding  the  efforts  required  for  model  solutions.  In  this  section, 
sensitivity  analysis  methods  for  the  space-time  problems  are  reviewed.  A 
computationally  efficient  method  in  calculating  sensitivity  coefficients, 
which  utilizes  the  information  obtained  in  solution  procedures,  is  also 
discussed. 


2.6.1  Description  of  Sensitivity  Coefficients 

For  a  space-time  problem  in  which  both  the  system  variables  and  the 
system  parameters  are  functions  of  space  and  time,  the  functional  derivative 
defined  below  provides  elementary  sensitivity  information 


6Ym  (x>t) 
(*' , t ' ) 


(2-55) 


The  sensitivity  density  represents  the  response  of  the  dependent 

variable  m  at  position  x,  and  time  t  with  respect  to  a  variation  of  the 
parameter  <t>  at  position  x'  and  time  t'.  Thus,  the  sensitivity  density  provide 
a  detail  system  response  to  the  variation  of  the  particular  parameter  value  at 
(x' ,  t').  However,  the  sensitivity  equations  should  be  solved  repeatedly  for 
each  variation  of  parameters.  Thus,  the  complete  sensitivity  information  for 
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a  single  parameter  may  require  solution  of  up  to  "n  x  m"  equations  (where  n 
is  the  number  of  nodes,  m  is  the  number  of  time  steps  taken).  Then,  the 
calculated  sensitivity  densities  may  be  examined  individually  or  be  integrated 
to  provide  a  total  variation  of  dependent  variable  of  interest  (6Ym), 

6Ym(x,  t;  xx  <  x  <  x2,  tt  <  t  <  t2)  = 


t2  X2 


J  J  fim  l  (X»ti  X'»  *')  fi^m(X'»  *'>  dX'  dt' 


fcl  X1 


(2-56) 


where  x^  and  X£  define  the  boundary  of  the  space  the  spatial  region,  and  t^ 
and  t2  define  the  temporal  duration. 

For  computational  efficiency,  the  sensitivity  coefficients  based  on 
constant  relative  variation  of  parameters  rather  than  the  sensitivity 
densities  are  used  in  this  study.  In  order  to  take  advantage  of  the  constant 
relative  variation  assumption,  the  parameter  (0£)  is  redefined  by  performing 
linear  transformation, 

*t  =  u  +  $£)*£  (2-57) 


where  is  the  constant  parameter  perturbed  in  the  sensitivity  analysis  and 
is  set  equal  to  zero  for  each  parameter.  Then  the  familiar  sensitivity 
coefficients  are  derived  by  taking  partial  derivative  with  respect  to 


3Y 


(2-58) 
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The  above  sensitivity  coefficients  may  be  divided  by  *m  to  provide  normalized 
sensitivity  information. 

It  is  often  of  interest  to  study  the  effects  of  variation  of  parameters 
in  finite  region  of  physical  coordinates  (i.e.,  x^  <  x  <  X2)  or  a  certain 
model  condition.  For  example,  many  of  thermal  properties  as  well  as  chemical 
kinetic  properties  at  high  temperatures  are  not  well  known.  Thus,  it  would  be 
helpful  to  study  the  importance  of  high  temperature  data  for  these  parameters. 
This  kind  of  sensitivity  information  can  be  obtained  by  applying  the  linear 
transform  defined  in  Eq.  (2-57)  only  to  the  region  of  interest.  If  the 
constant  variation  (not  constant  relative  variation)  is  assumed,  then  the 
following  transformation  is  used 

+  C  •  (2-59) 


The  sensitivity  coefficients  can  be  calculated  by  the  same  procedure  discussed 
above . 


2-6.2  Calculation  of  SensitivUy._c<?^ffi<?isnts 

The  sensitivity  coefficients  can  be  calculated  by  coupling  the  model 
equations  with  the  sensitivity  equations.  However,  this  approach  tends  to 
increase  the  numerical  stiffness  and  the  resulting  matrix  size  becomes  large 
especially  when  a  large  number  of  parameters  are  considered.  Here,  the  direct 
decoupled  method  (DDM)  proposed  by  Dunker^-^O  was  employed,  in  which  the 
sensitivity  equations  are  decoupled  from  the  model  equations,  and  integrated 
one  time  step  behind  the  integration  of  the  model  equations.  In  this  way,  the 
sensitivity  coefficients  are  efficiently  calculated  by  reusing  information 
obtained  from  solution  procedures. 
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Here,  the  DDH  is  extended  to  treat  the  implicit  form  of  ODEs  by  following 
Dunker's  work.  A  set  of  implicit  ODEs  rather  than  explicit  ODEs  are  resulted 
from  discretization  of  PDEs  by  using  the  moving  node  FEM,  which  can  be  written 
as 


A 


du 

dt 


=  g(u,  t,  <(> ) 


(2-60) 


The  Gear  algorithm  further  discretizes  the  above  ODEs  by  using  finite 
difference  approximations  to  yield 

i=q 

A(tn,u)  (u  -  £  °iyn-i^  "  h  &0  8^’“*  =  0  (2-61) 

i=l 

Where  q  is  the  order  of  integration  (here,  1  <  q  <  5)  and  0O  is  the 
coefficient  used  in  the  Gear  algorithm.  This  finite  difference  formulation 
may  be  re-written  in  a  more  convenient  form  for  a  modified  Newton  interaction 
technique. 


A(t  ,u)  (u  -  u“)  -  h  0  z(t  , 
n  non 


h,  u,  u  )  =  0 
*  n 


(2-62) 


Z(u)  =  g(t  ,u)  -  A(t  ,u) 
n  m 


u  -  a 

_U _ a 

hPo 


(2-63) 
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The  solutions  are  obtained  by  the  modified  Newton  iteration. 


[T<Vy>  '  h  6o  J(tn-y)l  (yi,m+1  '  y”>  ‘  ' 


h  3  z(t  , 
o  n’ 


m. 

yn> 
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The  sensitivity  equations  can  be  derived  by  applying  the  partial  derivative  to 
Eq.  2-10. 


i=q 

[A(tn,y)  -  h  3q  J(tn,y)]  Sn  =  A(tn,y)  £ 

i=l 

The  similarity  between  Newton  iteration  formula  for  solutions  and  sensitivity 
equations  become  quite  apparent  by  comparing  these  two  equations.  The  DDM 
greatly  increases  computational  efficiency  by  taking  advantage  of  the 
similarity  found  between  these  two  equations.  The  procedure  of  the  DDM  is  as 
follows:  the  sensitivity  equations,  Eq.  (2-65)  and  the  model  equations, 

Eq.  (2-64),  are  solved  sequentially.  The  model  equations  are  integrated  first 
from  tn_j  to  tn  by  using  the  Newton  iterative  formula.  During  this  solution 
procedure,  the  parameters  related  to  integration  such  as  a  time  step  (tn  - 
tn_i)  and  the  order  of  integration  (qn)  are  determined.  Then  the  Jacobian 
matrix  is  re-evaluated  and  LU  decomposed,  and  these  integration  parameters  as 
well  as  the  LU  decomposed  matrix  and  the  model  solutions  are  used  to  solve  the 
sensitivity  equations.  Each  sensitivity  equation  is  solved  separately.  This 
procedure  is  performed  iteratively  with  time. 
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3.0  HETEROGENEOUS  SURFACE  REACTION  MECHANISMS 


This  section  describes  the  heterogeneous  surface  reaction  mechanisms  that 
will  be  used  to  describe  particulate  boron  ignition,  oxidation,  and 
volatization.  Mechanisms  are  presented  for  two  idealized  cases  corresponding 
to  the  two  distinct  stages  of  boron  combustion  that  have  been  observed 
experimentally. 3-1, 3-2  The  first  reaction  set  includes  reactions  between  gas 
phase  species  and  a  boron  oxide  surface.  This  case  is  applicable  when  the 
boron  particle  has  a  significant  oxide  coating  (>  10  A)  and  is  important  in 
describing  the  chemical  transformations  that  promote  oxide  gasification.  The 
cond  case  treats  a  "clean"  boron  surface  (no  appreciable  oxide  coating)  is 
applicable  to  the  second-stage  ignition  and  full-fledged  boron  surface  burning 
that  occurs  at  high  temperatures  (ca.  1900  K)  once  the  oxide  layer  has  been 
removed. 

3. 1  Overview 

The  systematic  development  of  a  complete  and  self-consistent  kinetic 
model  for  the  heterogeneous  reactions  which  may  be  important  processes  in  the 
two  observed  stages  of  boron  combustion  proceeded  in  four  steps. 

First,  thermodynamic  equilibrium  calculations  for  several  hydrocarbon/ 
boron/air  mixtures  and  kinetic  calculations  using  the  homogeneous  gas  phase 
oxidation  model  developed  previous ly^”3 , 3-4  were  used  to  identify  potential 
gas  phase  reactants.  These  gas  phase  species  will  be  designated  as  Z(g). 

Second,  a  complete  heterogeneous  reaction  set  for  particulate  boron  and 
liquid  B2O3  gasification  in  a  hydrocarbon  combustion  environment  was 
formulated  in  terms  of  global  reactions  having  the  general  form 
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®2°3^  +  Zi/g^  +  Z2^8^  +  •••  +  zn(8) - >  gas  phase  products 


(3-1) 


and 


B(s)  +  ZL(g)  +  Z^(g)  +  ....  +  gn(g)  - >  gas  phase  products  .  (3-2) 


Third,  kinetic  mechanisms  consisting  of  surface  adsorption  and  desorption 
steps  were  developed  which  could  collectively  describe  the  global  reactions  in 
Eqs.  (3-1)  and  (3-2).  In  the  case  of  B2O3  gasification,  these  reaction  steps 
were  formulated  in  terms  of  the  structure  of  vitreous  boron  oxide  using 
elementary  thermochemical  and  kinetic  arguments.  In  the  case  of  particulate 
boron  surface  oxidation,  kinetic  mechanisms  and  rate  parameters  were 
formulated  on  the  basis  of  experimentally  observed  reaction  channels  and 
measured  reaction  cross  sections  for  small  boron  ion  clusters. 

Fourth,  kinetic  calculations  were  performed  and  model  results  analyzed 
for  internal  consistency  and  general  agreement  with  available  experimental 
data. 


Since  the  number  of  reactions  represented  by  Eqs.  (3-1)  and  (3-2)  is 
enormous,  it  was  necessary  to  selectively  prune  this  reaction  set  to  a  more 
computationally  manageable  number  of  reactions  without  excluding  important 
reaction  pathways.  Consequently,  thermochemical  and,  when  possible,  kinetic 
data,  in  conjunction  with  experimental  observations  and  earlier  modeling  work 
were  used  to  extract  from  Eqs.  (3-1)  and  (3-2)  the  most  important  reactions. 

For  a  first  order  representation  of  the  surface  kinetics,  three  criteria 
were  used  to  reduce  the  reaction  set  defined  by  Eqs.  (3-1)  and  (3-2).  First, 
only  reactions  whose  reaction  enthalpy  is  less  than  some  maximum  enthalpy  are 
included  in  the  reduced  reaction  set.  The  cutoff  enthalpy  was  defined  based 
on  the  vaporization  enthalpy  of  the  approriate  condensed  phase  to  ensure  that 
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the  surface  reactions  included  in  the  model  were  thermodynamically  competitive 
with  vaporization.  Second,  global  reactions  which  involve  gas  phase  reactants 
and/or  products  that  where  predicted  to  be  negligible  on  the  basis  of 
thermodynamic  and  kinetic  calculations  with  the  homogeneous  gas  phase 
oxidation  model  were  excluded  from  the  surface  reaction  model  to  ensure 
internal  consistency.  This  criterion  eliminated  global  reactions,  for 
example,  which  satisfied  the  first  criterion  but  yielded  B20(g)  as  a  final  gas 
phase  product.  Third,  only  reactions  which  are  first  order  in  gas  phase 
reactants  were  included  in  the  surface  kinetics  model  presented  here. 

This  criterion  was  initially  based  on  the  assumption  that  reaction  rates  for 
second  order  reactions  would  tend  to  be  slower  than  first  order  reaction  rates. 
This  assumption  was  found  to  be  consistent  with  model  calculations  (summarized 
in  Section  4)  for  high  temperature  boron  oxide  vaporization  and  boron 
combustion  which  indicated  that  only  ca.  1%  of  the  surface  was  covered  by 
adsorbed  species.  However,  the  reaction  mechanisms  adopted  in  this  model 
would  not  provide  an  adequate  description  of  oxide  layer  formation  at  lower 
temperatures  (T<  1600  K). 

It  is  important  to  note  that  these  criteria  are  not  rigid.  Rather,  they 
provide  defensible  reasons  for  the  particular  reactions  included  in  the 
kinetic  model,  yet  are  readily  expanded  to  include  additional  reactions  (by 
increasing  the  cut-off  enthalpy,  including  additional  species  in  the 
homogeneous  oxidation  model,  or  incorporating  second  order  processes)  as 
additional  experimental  or  theoretical  work  dictates  is  necessary. 

3.2  Gas  Phase  Boron  Speciation  and  Thermodynamics 

3.2.1  Gas  Phase  Reactants.  Heats  of  Formation  and  Entropies 

Thermochemical  equilibrium  calculations  and  results  from  the  application 
of  gradient  sensitivity  analysis  to  the  kinetic  calculations  for  the  homo¬ 
geneous  gas  phase  oxidation  indicate  that  H0B0(g),  6203(g)  and  B02(g)  ars 
major  gas  phase  species  in  oxygen  rich  environments  while  HBO(g),  8202(g)  and 
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BO(g)  are  the  primary  species  in  fuel  rich  environments. 3-3»3-4  Thus,  these 
species,  in  conjunction  with  several  hydrocarbon  combustion  intermediates 
and  products,  will  be  viewed  as  comprising  the  primary  set  of  gas  phase 
species  which  can  react  with  the  boron  or  boron  oxide  surface.  The  complete 
set  of  primary  gas  phase  species  and  their  entropies  and  heats  of  formation 
are  given  in  Table  1.  The  heat  of  formation  for  HBO(g)  is  based  on  the 
theoretical  calculations  of  Page.3"5  ^11  other  data  is  from  recent  JANAF 
listings . 


Table  1  -  AHf  £98  “d  S2gg  for  Primary  Gas  Phase  Speciation 
(±  uncertainties  in  parentheses) 


Species 

^®f ,298 
(kcal/mole) 

s298 

(cal/mole-deg) 

BO 

0.0 

(2.0) 

48.63 

(0.01)  ; 

bo2 

-  68.0 

(2.0) 

54.93 

(0.05) 

b2o2 

-109.0 

(2.0) 

57.98 

b2o3 

-199.8 

(1.0) 

[67.8 

(1.0)] 

HBO 

-  60.0 

(0.7) 

48.40 

(0.1) 

HOBO 

-134.0 

(1.0) 

[57.30] 

H 

52.1 

(0.001) 

27.39 

(0.004) 

0 

59.6 

(0.024) 

38.47 

(0.005) 

OH 

9.3 

(0.29) 

43.89 

(0.01) 

H20 

-  57.8 

(0.01) 

45.11 

(0.01) 

CO 

-  26.4 

(0.04) 

47.21 

(0.01) 

C02 

-  94.1 

(0.011) 

51.07 

(0.03) 
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3.2.2  Gas  Phase  Bond  Energies 


The  heats  of  formation  in  Table  1  were  used  to  determine  the  bond 
energies  given  in  Table  2.  These  bond  energies,  along  with  an  assumed  bond 
order,  were  subsequently  used  to  estimate  reaction  enthalpies  and  activation 
energies  for  the  elementary  reactions  defined  in  Section  3.3.  Note  that  not 
all  the  species  listed  in  Table  2  are  primary  gas  phase  species.  However, 
they  were  all  considered  in  estimating  the  reaction  enthalpies  of  elementary 
reactions  and  heats  of  formation  for  surface  speciation  and  are  listed  here 
for  completeness. 

3.3  Definition  of  Elementary  Reactions 

As  noted  in  Section  3.1,  kinetic  mechanisms  for  the  global  reactions 
represented  by  Eqs.  (3-1)  and  (3-2)  were  formulated  in  terms  of  simple 
adsorption  and  desorption  reaction  steps.  In  addition,  a  few  second  order 
reactions  were  adopted  to  ensure  that  the  global  reactions  were  reversible. 

As  a  result,  four  classes  of  reactions  were  defined  in  the  present  model  to  be 
elementary  and  have  reaction  rates  which  could  be  expressed  as  a  simple 
Arrhenius  function  of  temperature.  These  reactions  may  be  written,  in 
general,  as: 


S  +  Z(g)  <--->SZ(c)  (3-3) 

SZ(c)  - >  gas  phase  products  (3-4) 

SZ1(c)  +  SZ^fc)  — >  gas  phase  products  (3-5) 

SZ^c)  +  Z^(g)  — >  gas  phase  products  (3-6) 


where  S  is  a  reactive  surface  site  and  SZ(c)  is  a  surface  complex. 
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Table  2  -  Gas  Phase  Bond  Energies  (±  uncertainties  in  parentheses) 


Bond  Order 


(kcal/mole) 


HB=0 

B=0 

0BB=0 

OBOB-O 

HOB-O 

HO-BO 

(OH)2B-OH 

OBO-BO 

OB-O 

BO-B 

OB-BO 

(HO)2B-B(OH)2 

H-B 

HB-H 

H-BO 

H-H 

0=0 

H-0 

HO-H 

H-OBO 

C=0 

0C=0 


225.4 

(3.0) 

192.8 

(3.5) 

149.8 

(27.0) 

145.8 

(2.0) 

143.3 

(3.3) 

132.5 

(4.8) 

131.0 

(5.0) 

127.0 

(4.5) 

120.0 

109.0 

(6.0) 

76.0 

(13.0) 

78.2 

(0.04) 

110.3 

98.7 

(15.5) 

103.0 

(0.5) 

118.0  (0.1) 

101.0 

(0.3) 

118.0 

(0.1) 

118.1 

(3.0) 

256.0 

(0.04) 

127.6 

(0.06) 

4.28 

4.47 

5.12 

4.38 

5.12 

5.00 

11.11 

5.54 
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3.3.1  Adsorption  Reactions 


For  our  present  purposes,  adsorption  will  refer  strictly  to  chemisorption 
processes.  Ideally,  adsorption  channels  and  rate  parameters  for  the  kinetic 
model  would  be  based  on  experimental  data  or  microscopic  theoretical  calcula¬ 
tions.  Unfortunately,  this  approach  was  only  partially  feasible  in  the 
current  model  because  of  the  lack  of  experimental  data.  It  was  therefore 
necessary  to  formulate  adsorption  reactions  and  estimate  rate  parameters  using 
simple  thermodynamic  and  kinetic  arguments.  This  was  particularly  true  for 
the  B2O3  surface.  When  experimental  data  was  lacking,  therefore,  chemi¬ 
sorption  at  the  surface  of  a  condensed  phase  was  treated  as  follows: 


a.  The  clean  surface  was  assumed  to  have  a  particular  structure  in 
terms  of  which  it  was  possible  to  define  reactive  surface  sites  or 
reactive  surface  structural  units. 

b.  Potential  adsorption  channels  and  surface  speciation  was  defined  by 
considering  the  various  bonding  possibilities  between  a  surface  site 
and  a  gas  phase  species  Z. 

c.  Gas  phase  bond  energies  and  assumed  surface  species  structure  were 
used  to  estimate  adsorption  enthalpies.  The  most  thermodynamically 
stable  surface  complexes  that  were  consistent  with  the  overall 
global  reactions  in  Eqs.  (3-1)  and  (3-2)  were  adopted  as  the 
dominant  adsorption  channel. 

d.  The  adsorption  rate  constant  was  taken  to  have  the  form-^"? 


k  (T)  =  k  (T)  x  s  x  e(~Ea/kBT) 
a  max  o 

where 


(3-7a) 


k 


max 


[kBT/2irm]1/2 


(3.7b) 
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In  Eq.  (3.7)  is  the  number  of  molecules  of  mass  m  striking  the 

surface  per  unit  area  per  second,  sQ  is  the  sticking  probability, 

Ea  is  the  activation  energy  and  kg  is  Boltzmann's  constant.  If 
there  is  no  activation  energy  and  if  every  molecule  striking  the 
surface  chemisorbs  (sQ  =  1),  then  ka  =  k^ny  is  an  upper  limit  to 
rate  constant. 

e.  Lastly,  the  chemisorption  rate  is  given  by  the  Langmuir  rate^-? 

r  (T)  =  k  (T)  x  (1-0)  x  [Z(g) ]  (3.8) 

a  a 

for  immobile  adsorption,  and 

r  (T)  =  k  (T)  x  ( 1—0)2  x  (Z(g) ]  (3.9) 

a  a 

for  mobile  adsorption  where  8  is  the  fraction  of  surface  sites 
covered  by  absorbate. 

3.3.2  First  Order  Desorption  Reactions 

The  global  reactions  defined  by  Eqs.  (3-1)  and  (3-2),  may  be 
characterized  by  an  initial  adsorption  of  gas  phase  Z(g)  to  yield  the  surface 
complex  SZ(c),  followed  by  decomposition  of  this  complex  to  yield  the  final 
gas  phase  product  Z'(g).  A  complete  reaction  sequence  can  be  written  as 

k  (T)  k  (T) 

sn  +  z(s>  <!:::>  snz(c)  “ — >  sn-i  +  z'(g)  (3”10) 


The  desorption  ethalpy  is  then  given  by 


H  =  H  -  H 
d  r  a 


(3-11) 
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where  Ha  is  the  estimated  adsorption  enthalpy  and  Hr  is  the  reaction  enthalpy 
determined  from  heats  of  formation  for  reactants  and  products.  The 
first-order  desorption  rate  constant  (probability  per  sec)  is  given  by 


kd(T)  =  Ad 


(-Ed/kBT) 

x  e 


(3-12) 


As  a  first-order  approximation,  the  coefficient  A^  was  expressed  as  the 
transition  state  frequency  factor  kgT/h  and  the  activation  energy  was 
estimated  using  gas  phase  bond  energies. 

The  desorption  reaction  rate  for  the  inverse  adsorption  step,  i.e., 

SnZ(c)  - >  Sn  +  Z(g)  (3-13) 

was  estimated  from  the  adsorption  rate  constant  ka  and  equilibrium  constant 
Keq-  Keq  was  computed  using  the  estimated  adsorption  enthalpy  and  tabulated 
values  for  gas  phase  entropies  from  Table  1. 

3.3.3  Second  Order  Reactions 

For  the  B2O3  surface,  the  simplest  adsorption  model  was  typically  found 
to  result  in  the  formation  of  two  surface  complexes  for  each  gas  phase 
reactant  adsorbed.  Consequently,  for  the  global  reactions  represented  by  Eq. 
(3-1)  to  be  reversible  it  was  necessary  to  include  the  second  order  desorption 
steps  represented  by 

SZ^(c)  +  SZ^fc)  - >  gas  phase  products  (3-14) 
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in  the  kinetic  model.  Similarly,  it  was  found  that  for  the  global  reactions 
represented  by  Eq.  (3-2)  for  particulate  boron  to  be  reversible,  the  kinetic 
model  had  to  include  second  order  reactions  of  the  form 

SZj(c)  +  Z^(g)  — >  S  +  gas  phase  products  .  (3-15) 


In  both  these  cases,  rate  parameters  were  estimated  from  Keq  for  the  global 
reaction  and  the  adsorption/desorption  rates  defined  in  Sections  3.3.1  and 
3.2.2. 


3.4  Boron  Oxide  Surface  Reactions 

Oxidation  reactions  between  oxygen  and  a  clean  boron  surface  are 
typically  exothermic  (see  Section  3.5.1).  Subsequent  reactions,  however, 
between  gas  phase  reactants  and  the  oxide  layer  that  forms  on  the  boron  sur¬ 
face  are  typically  highly  endothermic.  Consequently,  oxide  layer  formation 
tends  to  inhibit  oxidation  until  the  particle  temperature,  through  convective 
and  radiative  heat  transfer  or  self-heating,  reaches  the  vaporization 
temperature  of  B2O3.  Removal  of  the  oxide  layer  can  be  aided  by  surface 
reactions  which,  although  endothermic,  require  less  energy  than  is  needed  for 
vaporization.  Experimental  evidence  indicates  that  hydrogen  containing 
species  tend  to  increase  B2O3  vaporization  rates-*"®  and  to  decrease  both  the 
ignition  temperature  and  ignition  time  delay  for  boron  particles. 3-1, 3-2 

3.4.1  Global  Reactions 

Table  3  lists  reactions  between  condensed  B2O3  and  primary  gas  phase 
species  to  produce  gas  phase  products.  Following  the  criteria  discussed  in 
Section  3.1,  only  reactions  which  are  first  order  in  gas  phase  reactants  and 
whose  reaction  enthalpy  is  less  than  110  kcal/mole  are  included;  reactions 
with  a  larger  reaction  enthalpy  would  be  competitive  with  B2O3  vaporization 
(AH  -100  kcal/  mole).  Also,  reactions  which  involved  B(0H)2(g)  as  a  final 
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product  were  excluded,  even  if  the  reaction  enthalpy  was  less  than 
110  kcal/mole,  because  equilibrium  calculations  for  the  oxidation  of  a 
JP4/solid  boron  mixture  in  air  indicated  that  B(0H)2  would  be  a  minor 
species. 3-3, 3-4  ^s  shown  in  Table  3,  there  are  only  three  reactions  which. 


Table  3  -  Boron  Oxide  Surface:  Global  Reactions  (Excludes  Reactions 
With  Reaction  Enthalpies  >  110  kcal/mole  and  Reactions 
Involving  B(0H)2(g)  As  a  Product) 


Reaction 

AH 

(kcal/mole) 

B203( 1)  +  0(g)  - >  2B02(g) 

B203(1)  +  0H(g)  - >  H0B0(g)  +  B02(g) 

B203(1)  +  H20(g)  - >  2H0B0(g) 

108.4  (4.5) 

92.7  (3.8) 

93.8  (2.5) 

based  on  current  thermochemical  data,  satisfy  these  constraints.  The 
important  gas  phase  reactants  are  the  hydrocarbon  combustion  products  0(g), 
0H(g)  and  H20(g),  although  0(g)  is  not  expected  to  have  much  influence  on  the 
gasification  process. 

3.4.2  Oxide  Structure 

Structure  of  vitreous  boron  trioxide  has  been  a  subject  of  considerable 
disagreement  and  controversy.  Early  concepts  ranged  from  that  of  a  molecular 
liquid  consisting  of  the  dimers,  B^Og,  an  ionic  melt  involving  complex 
boron-oxygen  ions,  to  that  of  a  highly  associated  network  of  interlinking  BO3 
triangles . ^“9  The  latter  structure  is  in  general  agreement  with  the  classical 
theory  of  the  vitreous  state  for  B2O3,  Si02  and  Ge02- 

Early  ^B  NMR  spectra  were  thought  to  generally  support  the  model  of  a 
random  network  of  planar  BO3  triangular  units  interconnected  through  the 
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oxygen  atoms  at  the  vertices . 3-10,  3-11  However,  a  double  bond  resonance  in 
which  boron  attracts  one  electron  from  oxygen  and  forms  a  predominately  ionic 
bond  in  addition  to  the  covalent  sigma  bond  was  introduced  to  bring 
experimental  measurements  of  the  quadrupole  coupling  constant  into  agreement 
with  trigonal  BO3  structural  units. 

The  most  recent  studies  based  on  infrared  and  Raman  spectroscopy, 3-12 
X-ray, 3-13  NMR  techniques 3” 14  have  been  interpreted  to  support  the 
existence  of  a  random  network  of  boroxol  rings  interconnected  by  a  random 
number  of  BO3  triangles.  NMR  *?0  spectra  for  B2O3  glass,  for  example, 
indicated  the  presence  of  two  distinct  oxygen  sites  corresponding  to  oxygen 
atoms  inside  and  outside  the  boroxol  ring. *4  Spectra  of  oxygen  atoms 
outside  the  boroxol  rings  were  characterized  by  much  larger  distributions  in 
the  coupling  constant  and  asymmetry  parameter  than  the  spectra  for  oxygen 
atoms  in  the  boroxol  rings.  In  addition,  an  X-ray  diffraction  study  of  molten 
boron  trioxide  (650  C)  concluded  that  the  melt  contained  a  mixture  of  boroxol 
rings  linked  by  independent  BO3  units. 3-15 

Based  on  these  latter  experimental  results,  the  structure  of  condensed 
phase,  vitreous  B2O3  can  be  represented  as  shown  in  Figure  3.1a.  However,  for 
the  purpose  of  estimating  adsorption  and  desorption  kinetic  and  thermochemical 
parameters  for  the  global  surface  reactions  listed  above,  we  can  initially 
neglect  subtle  differences  in  the  bond  angles  and  bond  strengths  of  different 
surface  sites,  since  the  uncertainty  in  the  estimated  reaction  rates  should 
encompass  any  differences  in  surface  site  reactivity  due  to  the  presence  of 
broxol  rings.  Thus,  for  the  present  model,  a  reactive  surface  site  will  be 
represented  as  shown  in  Figure  3.1b.  Our  model  for  condensed  B2O3  is, 
consequently,  more  consistent  with  studies  based  on  computer  modeling^-***  and 
molecular  dynamics  calculations^-*?  which  hold  that  the  structure  of  vitreous 
boron  trioxide  can  be  understood  on  the  basis  of  a  continuous  random  network 
of  BO3  triangles  containing  no  ring  structures. 
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Figure  3.1.  (a)  One  proposed  structure  of  vitreous  B2O3  consisting  of 

boroxyl  rings  interconnected  by  BO3  units,  (b)  Surface 
structural  unit  used  in  the  kinetic  model  to  describe 
surface  adsorption  and  desorption  reaction  steps. 
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3.4.3  Adsorption  Thermodynamics  and  Kinetics 


For  the  B2O3  surface,  the  basic  structural  unit  used  to  define  surface 
spec iat ion  is 

°>B-0-B<0 

where  the  oxygen  atoms  are  bonded  through  B-0  bonds  in  trigonal  BO3  units. 
Surface  speciation  is  then  defined  by  treating  the  reactions  between  gas  phase 
reactants  and  the  B2O3  reactant  as  a  simple  'bimolecular '  reaction.  Based  on 
the  reactions  in  Table  3,  0(g),  0H(g)  and  H20(g)  are  the  only  gas  phase 
reactants  which  need  to  be  considered  at  this  stage. 

For  each  of  the  three  gas  phase  reactants  (0,  OH,  and  H2O)  in  which  we 
are  primarily  interested,  there  are  several  potential  surface  complexes  that 
could  be  formed  upon  adsorption.  However,  it  was  found  that  most  of  these 
complexes  were  either  more  than  10  Kcal  endothermic,  whereas  the  most 
thermodynamically  favored  complexes  were  typically  1-20  kcal  exothermic,  or 
would  require  significant  bond  rearrangement  to  yield  final  gas  phase  products 
consistent  with  the  global  reaction  products  in  Table  3.  Since  experimental 
data  for  the  kinetic  mechanisms  is  not  currently  available,  the  most 
straightforward  procedure  was  to  adopt  only  the  most  thermodynamically  stable 
surface  complexes,  based  on  estimated  adsorption  enthalpies,  that  were 
consistent  with  the  global  surface  reactions  in  Table  3.  As  a  result,  the 
following  three  adsorption  reactions  were  used  in  the  kinetic  model: 

°>B-0-B<0  +  0(g)  — >  °>B-0*  +  °>B-0  (3-16) 

°>B-0-B<0  +  0H(g)  — >  °>B-0’  +  °>B-0H  (3-17) 

°>B-0-B<0  +  H20(g)  — >  °>B-0H  +  °>B-0H  (3-18) 


3-14 


In  Eqs.  3-16  and  3-17,  0*  represents  oxygen  with  an  unpaired  electron. 

Each  of  the  adsorption  reactions  can  be  viewed  as  resulting  from  the 
fragmentation  of  a  >B-0-  bond  on  the  surface  and  the  formation  of  a  new  bond 
between  the  electron  deficient  boron  and  the  more  electronegative  portion  of 
the  gas  phase  reactant. 

Estimated  chemisorption  enthalpies  and  rate  parameters  for  the  adsorption 
reactions  are  summarized  in  Table  4.  The  adsorption  enthalpies  and  activation 
energies  are  based  on  gas  phase  bond  energies.  The  sticking  probabilities  are 
based  on  a  rough  analogy  to  the  adsorption  probabilities  for  these  reactants 
with  other  substrates. 


Table  4  -  Boron  Oxide  Surface:  Adsorption  Reactions 


Z(g) 

AH29b 

(kcal/mole) 

Min  Est  Max 

MT)  = 

^max 

(cm/s) 

km ax  x  s0  x  e(— Ea/RT) 

so  Ea 

(kcal/mole) 

Adsorption 

Products 

0 

0 

-10 

-50 

925xT0-5 

0.2  -  0.7  0  -  30 

2 [>B0. ] 

OH 

0 

-10 

-25 

897xT0-5 

0.2  -  0.7  0  -  30 

[>B0H]  +  [>B0.] 

h2o 

0 

-10 

-35 

872xT0-5 

0.1  -  0.6  5  -  50 

2 [>B-0H] 

3.4.4  Desorption  Reactions 

The  surface  species  that  result  from  the  adsorption  reactions, 

Eqs.  (3-16)  and  (3-17),  may  be  written  as  >B-X  with  X=(0H,0.)  where  the  dot 
indicates  a  lone  electron.  Then  gas  phase  products  consistent  with  the 
reactions  in  Eq.  (3.1)  can  be  viewed  as  arising  from  the  following  sequence  of 
bond  cleavage  and  formation: 
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0 

0 


/B°2 


\ 


B-X 


/ 

\B0 


- >  bo2-o.  +  .bo2 


J' 


)B-X  - >°>B-0-B<°  +  0=B-X 

(3-19) 


It  is  clear  that  this  reaction  sequence  may  not  be  elementary. 
Nevertheless,  it  will  be  assumed  that  reaction  of  the  two  radicals  to  yield 
the  B2O3  surface  structural  unit  and  formation  of  the  0=B-  double  bond  is  much 
faster  than  fragmentation  of  the  two  B-0  bonds.  The  decomposition  reaction 
can  then  be  written  as 


(bo3)2-b-x 


0=B-X 


(3-20) 


In  Equation  (3-20)  k^  is  the  first-order  desorption  rate  constant  given  by 


=  Al 


(3-21) 


where  the  desorption  coefficient  is  given  approximately  by  the  transition 
state  frequency  factor  kgT/h.  Desorption  reaction  products,  reaction 
enthalpies  and  reaction  rates  are  summarized  in  Table  5. 


Table  5  -  Boron  Oxide  Surface:  First  Order  Desorption  Reactions 


Reaction 


h298 

(Kcal/mole) 

k=A1*e(-E1/RT) 

Min  Est  Max 

Al, 

(ps"1) 

El 

(Kcal/mole) 

>B-0H  ->  H0B0(g) 
>B0.  ->  B02(g) 


39  51  63 
42  54  66 


20 


40  -  60 


42 


54 
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20 


40  -  70 


It  can  be  seen  that  these  reactions  in  conjunction  with  the  chemisorption 
reactions  in  Table  4  can  be  used  to  describe  the  reactions  in  Eq.  (3-1)  for 
the  B2O3  surface  which  are  first  order  in  gas  phase  reactants.  For  example, 
the  reaction  between  the  B2O3  surface  and  the  hydroxyl  radical  OH(g)  to  yield 
HOBO(g)  and  B02(g)  can  be  describe  by  the  following  reaction  sequence: 


°>B-0-B<q  +  OH(g)  - >°>B-OH  +°>B0. 

(3-22) 

°>B-OH  <=-=>  HOBO(g) 

(3-23) 

/ — V 
60 
w 
CM 
O 
OQ 

A 

1  1 

1  1 

1  1 

V 

d 

oa 

A 

O 

(3-24) 

3.5  Boron  Particle  Surface  Reactions 

Once  the  oxide  layer  is  removed,  a  boron  particle  burns  rapidly  with  a 
surrounding  flame  sheet  which  is  detached,  but  near  the  surface  (see 
Figures  1.1  and  1.2).  Since  the  volatility  of  boron  is  low,  boron  burns  more 
like  a  carbon  particle  rather  than  a  liquid  hydrocarbon  droplet. 3~18  jn 
particular,  gasification  of  the  boron  particle  occurs  by  chemical 
transformation  into  more  volatile  components.  Thus,  heterogeneous  and 
homogeneous  chemistry  play  a  more  dominant  role  in  the  gasification  of  boron 
than  for  liquid  hydrocarbons.  For  example,  liquid  hydrocarbon  droplets  gasify 
primarily  by  a  simple  vaporization  process  which  proceeds  in  the  presence  or 
absence  of  an  oxidant. 

The  volatile  component  formed  during  boron  gasification  was  initially 
believed  to  be  BO(g),  however  there  may  be  significant  levels  of  HOBO(g)  and 
HBO(g),  as  well  as  8203(g),  8202(g)  and  B02(g).  It  is  interesting  to  note 
that  the  dominant  equilibrium  gas  phase  boron  species  HOBO(g)  and  HBO(g)  (in 
boron/JP4/air  systems)  are  isoelectronic  with  the  primary  hydrocarbon 
oxidation  products  C02(g)  and  CO(g),  respectively. The  formation  of 
BO  and  other  suboxides  and  suboxyhydrides  may  be  compared  to  the  formation  of 
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CO  during  carbon  gasification.  There,  C02(g),  the  major  carbon  combustion 
product,  diffuses  back  to  the  carbon  particle  surface  to  react  heterogeneously 
to  form  more  CO(g).  Boron  oxides  and  oxyhydrides  may  play  similar  roles  in 
boron  combustion. 

3.5.1  Global  Reactions 

Table  6  lists  reactions  between  solid  boron  and  primary  gas  phase  species 
to  produce  gas  phase  products.  Only  reactions  whose  reaction  enthalpy  was 
less  than  150  kcal/mole  are  included;  more  endothermic  reactions  would  not  be 
thermodynamically  competitive  with  boron  vaporization  since  the  heat  of  forma¬ 
tion  of  B(s)  is  approximately  133  kcal/mole.  In  addition,  this  reaction  list 
is  limited  to  those  reactions  which  are  first  order  in  gas  phase  reactants. 
Lastly,  reactions  which  involved  B20(g),  boron  hydrides  and  species  such  as 
B(0H)2(g),  H3B03(g)  and  1138303^)  were  excluded  because  these  species  are 
minor  gas  phase  constituents . 3“3>  3”^ 

3.5.2  Experimental  Rate  Data 

Two  sets  of  experimental  data,  summarized  below,  were  used  to  estimate 
reaction  steps  and  rate  parameters  for  the  global  reactions  in  Table  6. 

Cluster  Ion  Reactions 

Anderson  et.  al3-19-3-22  have  performed  a  series  of  investigations  into 
the  reactivity  of  boron  ion  clusters  using  the  cluster  ion  beam  technique.  In 
this  work,  the  reactions  of  thermalized,  size  selected  Bj_i3+  ions  were 
examined  under  single  collision  conditions  for  several  gas  phase  reactants, 
including  02(g),  D20(g),  C02(g)  and  C0(g).  Using  experimental  reaction  cross 
sections,  temperature  dependent  rate  constants  for  these  reactants  were  then 
obtained  by  dividing  the  measured  cross  section  by  a  collision  cross  section 
and  averaging  over  a  Boltzmann  distribution.  For  reference,  the  Bj2+  cluster 
reactions  and  product  channels  observed  by  Anderson's  group  are  listed  in 
Table  7. 
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Table  6  -  Boron  Surface:  Global  Reactions 


Reaction 

AH(298) 

(Kcal/mole) 

B(s)  +  0(g)  ->  B0(g) 

-59.6  (2.0) 

B(s)  +  02(g)  -*•  B02(g) 

-68.0  (2.0) 

B(s)  +  02(g)  ->  B0(g)  +  0(g) 

59.6  (2.0) 

2B(s)  +  02(g)  -»  8202(g) 

-109.0  (2.0) 

2B(s)  +  02(g)  +  2B0(g) 

0.0  (4.0) 

B(s)  +  0H(g)  -  HBO(g) 

-56.7  (1.0) 

B(s)  +  0H(g)  -  B0(g)  +  H(g) 

42.8  (2.3) 

B(s)  +  H20(g)  -  B0(g)  +  H2(g) 

57.8  (2.0) 

B(s)  +  H20(g)  HBO(g)  4-  H(g) 

62.5  (0.7) 

B(s)  +  B02(g)  -  B202(g) 

-41.0  (4.0) 

B(s)  +  B02(g)  -  2B0(g) 

68.0  (6.0) 

B(s)  +  B203(g)  ->  B202(g)  +  B0(g) 

90.8  (6.0) 

B(s)  +  H0B0(g)  -►  B0(g)  +  HBO(g) 

86.6  (3.7) 

B(s)  +  H0B0(g)  -  B202(g)  +  H(g) 

77.1  (3.0) 

B(s)  +  C02(g)  -  B0(g)  +  C0(g) 

67.7  (2.0) 

Of  the  15  reactions  in  Table  6  one  reaction  can  be  neglected  based  on 
Anderson's  boron  ion  cluster  work.  Specifically,  for  clusters  larger  than 
Bg+,  the  D2O  product  distributions  were  observed  to  be  dominated  by  a  single 
channel:  Bn+  +  D20(g)  -*•  Bn_]D+  +  DBO(g).  Accordingly,  we  can  assume  that 
HBO(g)  +  H(g)  is  the  primary  product  channel  for  the  reaction  between 
particulate  boron  and  H2O  and  neglect  the  channel  corresponding  to 
80(g)  +  H2(g) . 
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Table  7  -  Bi2+  Cluster  Reaction  Product  Channels 


Reaction 

Gas  Phase  Products 
(most  stable) 

Reaction  Probability 

1750  K 

2500  K 

b12+  +  02  Bn0+ 

BO 

0.00834 

0.008001 

B10+ 

b2°2 

0.046932 

0.052014 

-  Bll+ 

bo2 

5.3x10-9 

3.6x10-7 

b12+  +  d2°  b10d2+ 

b2o 

0.030898 

0.022496 

-  BllD+ 

DBO 

0.449861 

0.408516 

B^2+  C02  Bi20+ 

CO 

0.000321 

0.000528 

b12+  +  CO  ■+  B12C0+ 

— 

0.000137 

0.00016 

Note  also  that  the  reaction  probability  for 

B12+  +  02 - >  Bu+  +  B02(g)  (3-25) 

is  on  the  order  of  10-7  at  T=2500.  As  a  consequence,  in  the  present  model  we 
shall  assume  that  the  2nd  global  reaction  in  Table  6,  B(s)  +  02(g)  ■+•  B02(g), 
proceeds  via  the  reaction  sequence. 

B12+  +  02  — >  Bu0+  +  B0(g)  (3-26) 

Bn0+  +  B0(g) - >  B1q+  +  B02(g)  .  (3-27) 
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It  should  also  be  noted  that  in  the  channel  producing  8202(g)  no  surface 
complex  is  formed.  In  our  model,  however,  we  will  assume  that  for  larger 
particles  thermal  accommodation  will  allow  complex  formation  with  a  residence 
time  governed  by  the  desorption  rate  for  8202(c)  ■+  8202(g)*  However,  to 
maintain  internal  consistency  with  other  reactions,  additional  desorption 
channels  have  also  been  included. 

Lastly,  to  ensure  that  global  reaction  3  in  Table  8  is  reversible,  the 
reaction 


80(c)  +  0(g)  - >  B(s)  +  0£(g) 


(3-28) 


must  be  included  in  the  model.  Written  as  an  adsorption  reaction,  this 
incorporates  an  02(g)  reaction  channel  that  was  not  observed  by  Anderson. 

Flow  Reactor  Observations 

Rosner  et.  al. 3”24-3-26  have  investigated  the  heterogeneous  kinetics  for 
the  high  temperature  gasification  of  solid  boron  utilizing  low-pressure 
transverse  filament  flow  reactor  techniques  and  'real-time'  microwave  induced 
plasma  emission  spectroscopy  (MIPES)  element  detection.  This  work  included 
studies  of  reactions  between  solid  boron  and  0(g),  02(g),  C02(g),  H20(g)  and 
8203(g)*  Experimental  conditions  typically  covered  temperatures  between 
1300  K  and  2100  K  and  reactant  partial  pressures  between  ca.  10"^  and  10"1  Pa. 
Measured  reaction  rates  were  reported  in  terms  of  a  non-dimensional,  'overall' 
rate  constant  that  was  defined  as  the  ratio  between  the  net  flux  of  boron 
atoms,  irrespective  of  speciation,  emerging  from  the  surface  as  a  result  of 
chemical  reaction,  to  the  arrival  flux  of  gaseous  reactant.  The  reported 
reaction  probabilities  are  therefore  independent  of  reaction  products. 

In  the  present  model,  this  experimental  work  was  used  to  estimate  rate 
constants  for  boron  surface  reactions  with  0(g)  and  8203(g).  For  8203(g),  the 
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experimental  data  indicate  a  remarkably  high  reaction  probability  with  a 
maximum  close  to  unity  near  2000  K,  much  higher  than  that  for  the  02(g)/B(s) 
reaction  and  comparable  to  that  observed  for  atomic  reactants  such  as  0(g). 

Arrhenius  plots  indicate  that  the  reaction  is  not  elementary  over  a  broad 
range  of  conditions  and  that  there  may  be  several  different  rate-controlling 
reactions  with  rather  different  activation  energies.  Since  there  is  currently 
not  enough  experimental  data  to  speculate  on  the  reaction  mechanism,  we 
adopted  a  rather  crude  fit  to  the  experimental  data  which  represents  an 
approximate  average  reaction  probability  for  T  >  1800  K. 

3.5.3  Kinetic  Model 

A  kinetic  model  was  formulated  to  describe  the  global  reactions  given  in 
Table  6  in  terms  of  adsorption  and  desorption  reaction  steps.  Adsorption 
steps  were  generally  defined  by  analogy  to  observed  reaction  channels  for 
thermalir'.ed,  size  selected  boron  ion  clusters  under  single  collision 
conditions.  Adsorption  enthalpies  and,  in  some  cases,  desorption  activation 
energies  have  been  defined  in  terms  of  the  heats  of  formation  for  surface 
complexes.  These  were  estimated  using  gas  phase  bond  energies  as 


Hj  =  heat  of  formation 
H£  =  heat  of  formation 
H3  =  heat  of  formation 
H4  =  heat  of  formation 
H5  =  heat  of  formation 


of  B0(c)  --->  -71 

of  BH(c)  - >  26 

of  HB0(c)  — >  -45 

of  BCO(c)  — >  -80 

of  B£02(c)  — >  -70 


0  Kcal/mole  (3-29) 
0  Kcal/mole  (3-30) 
0  Kcal/mole  (3-31) 
0  Kcal/mole  (3-32) 
0  Kcal/mole  (3-33) 
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Adsorption  reactions  are  listed  in  Table  8.  Adsorption  enthalpies  are 
expressed  in  terms  of  estimated  heats  of  formation  for  the  appropriate  surface 
complex  and  tabulated  JANAF  heats  of  formation  for  the  other  reactants.  Rate 
parameters  are  based  on  Anderson's  Bj2+  reaction  cross  sections  and  Rosner's 
MIPES  studies.  All  these  reactions  are  reversible  and  rate  constants  for  the 
reverse  reaction  were  computed  in  the  kinetics  code  using  keq  and  ka. 

For  the  reactants  02(g),  H20(g),  C02(g)  and  CO(g)  we  have  assumed  that 
the  dominant  product  channels  correspond  to  those  observed  by  Anderson's  group 
for  Bj2+  clusters.  The  reaction  probability  per  collision  is  then  given  by 
the  measured  cross  section  divided  by  a  collision  cross  section  and  averaged 
over  a  Boltzmann  distribution.  The  temperature  dependent  reaction  probability 
is  then  fit  to  an  Arrhenius  function,  s=sQ  x  exp(-E/RT). 

Anderson's  work  could  not  determine  neutral  gas  phase  products. 
Consequently,  when  more  than  one  product  was  possible  [e.g.  2B0(g)  versus 
82^2(8)]  the  most  stable  products  were  used. 

The  last  three  reactions  in  Table  8  represent  adsorption  reactions  for 
gas  phase  reactants,  which  are  formed  by  surface  reactions.  As  in  the  model 
for  B2O3  surface  reactions,  they  are  included  to  ensure  reversibility  and  the 
rates  are  determined  from  Keq  for  the  appropriate  global  reactions. 

The  desorption  steps  needed  to  describe  the  global  reactions  given  in 
Table  6  are  listed  in  Table  9  along  with  estimates  of  the  desorption  enthalpy 
and  rate  parameters.  Since  these  desorption  reactions  are  all  first  order  in 
the  surface  complex,  the  desorption  rate  is  taken  to  have  the  form 

rd  =  Ad  x  T  x  exp(-Ed/RT)  .  (3-34) 
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Table  8  -  Boron  Surface:  Adsorption  Reactions 


Reaction 

km  ax 
(cm/sec) 

sO 

Ea 

(cal/mole) 

Ha 

(Kcal/mole) 

B(s)  +  H(g)  =  BH(c) 

3623 

0.1 

0.0 

H2-52. 1 

B(s)  +  0(g)  =  B0(c) 

909 

0.8 

0.0 

Hj-59.6 

2B(s)  +  02(g)  =  B202(c) 

643 

0.064 

1014.0 

-109.0 

2B(s)  +  02(g)  =  B0(c)  +  B0(g) 

643 

0.0082 

0.0 

«l 

B(s)  +  02(g)  =  B0(c)  +  0(g) 

643 

0.001 

1000.0 

Hj+59 . 6 

B(s)  +  0H(g)  =  B0(c)  +  H(g) 

882 

0.02 

1000.0 

Hj+42.8 

B(s)  +  0H(g)  =  HBO(c) 

882 

0.02 

1000.0 

H3-9.3 

2B(s)  +  H20(g)  =  BH(c)  +  HBO(g) 

857 

0.0062 

8000.0 

H2-2.2 

B(s)  +  B02(g)  =  B202(c) 

556 

0.005 

5000.0 

H5+68.O 

B(s)  +  B02(g)  =  B0(c)  +  B0(g) 

556 

0.005 

5000.0 

Hx+68.0 

B(s)  +  B203(g)  =  B202(g)  +  B0(g) 

436 

0.87 

0.0 

90.8 

B(s)  +  HOBO(g)  =  HBO(c)  +  B0(g) 

549 

0.0025 

4000.0 

H2+134. 0 

B(s)  +  HOBO(g)  =  B202(c)  +  H(g) 

549 

0.0025 

4000.0 

H5+I86.I 

B(s)  +  C02(g)  =  B0(c)  +  C0(g) 

548 

0.0013 

4517.0 

Hi+67.7 

B(s)  +  C0(g)  =  BCO(c) 

687 

0.00015  0.0 

H4+26.4 

B(s)  +  B0(g)  =  B20(c)  =  B0(c) 

Rate  computed 

from  Keq 

»1 

B(s)  +  B202(g)  =  B302(c)  =  B202(c) 

Rate  computed 

from  Keq 

H5+IO9.O 

B(s)  +  HBO(g)  =  HB20(c)  =  HBO(c) 

Rate  computed 

from  Keq 

H3+6O.O 
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The  desorption  prefactor  (A^xT)  was  taken  to  equal  the  transition  state 
frequency  factor  (kjjT/h)  and  the  activation  energy  as  set  equal  to  the 
endothermic  desorption  enthalpy.  Rate  parameters  are  not  given  for  those 
reactions  which  are  the  reverse  of  the  adsorption  reactions  listed  in  Table  10 
since  they  follow  from  directly  from  the  equilibrium  constant  and  are  computed 
during  execution  of  the  kinetics  code.  They  are  included  in  this  Table  only 
to  illustrate  all  the  desorption  channels  included  in  the  model. 


Table  9  -  Boron  Surface:  First  Order  Desorption  Reactions 


Reaction 

Hd 

(Kcal/mole) 

Ad 

(ps-1) 

Ed 

(Kcal/mole) 

BH(c)  -v  B(s)  +  H(g) 

52.1-H2 

Reverse  adsorption 

B0(c)  -  B0(g) 

-HI 

0.02 

-Hi 

B0(c)  ->  B(s)  +  0(g) 

59.6-H1 

Reverse  adsorption 

B202(c)  -*  B202(g) 

-109.0-H5 

0.02 

-109.0-H5 

-  2B(s)  +  02(g) 

-H5 

Reverse  adsorption 

-  B(s)  +  B02(g) 

-68.0-H5 

Reverse  adsorption 

-  2B0(g) 

-H5 

Neglect 

HBO(c)  -►  HBO(g) 

-60.0-H3 

0.02 

-60.0-H3 

-  B0(g)  +  H(g) 

52.1-H3 

0.02 

52.1-H3 

-*•  B(s)  +  0H(g) 

9 . 3-H3 

Reverse  adsorption 

BCO(c)  ->  B(s)  +  C0(g) 

-26.4-H4 

Reverse  adsorption 

There  are  potentially  at  least  four  major  channels  associated  with  the 
desorption  of  B202(c),  two  of  which  are  the  reverse  of  adsorption  reactions. 
The  remaining  two  channels  correspond  to  the  product  distributions  B202(g)  and 
2B0(g),  with  the  former  being  approximately  109  Kcal/mole  more  stable.  The 
thermodynamically  favoured  B202(g)  was  taken  as  the  dominant  channel. 
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To  ensure  that  the  global  reactions  in  Table  6  were  reversible,  a  few 
second  order  reactions  were  in  the  kinetic  model.  All  these  reactions  are 
reverse  of  adsorption  reaction  in  Table  8  so  that  in  model  calculations  the 
rate  constant  is  obtained  directly  from  Keq.  These  reactions  are  listed  in 
Table  10  to  help  clarify  specifically  which  second  order  reactions  are 
currently  included  in  the  model. 


Table  10  -  Boron  Surface:  Second  Order  Reactions 


Reaction 


B0(c)  +  H(g) 
B0(c)  +  0(g) 
B0(c)  +  C0(g) 
B0(c)  +  B0(g) 
B0(c)  +  B0(g) 
HBO(c)  +  B0(g) 
B202(c)  +  H(g) 
BH(c)  +  HBO(g) 


-►  B(s)  +  0H(g) 

-*•  B(s)  +  02(g) 

-  B(s)  +  C02(g) 

-  2B(s)  +  02(g) 

-  B(s)  +  B02(g) 
B(s)  +  H0B0(g) 

-♦  H0B0(g) 

-  2B(s)  +  H20(g) 


3 . 6  References  for  Section  3.0 


3-1  Macek,  A.,  "Combustion  of  Boron  Particles  at  Atmospheric  Pressure," 
Comb.  Sci.  and  Tech.  2,  181  (1969). 

3-2  Macek,  A.,  "Combustion  of  Boron  Particles:  Experiment  and  Theory," 

Fourteenth  Symp.  (Int.)  on  Comb.,  pg.  1401,  The  Combustion  Institute, 
Pittsburgh,  PA  (1972). 

3-3  Yetter,  R.A.,  Rabitz,  H. ,  Dryer,  F.L.,  Brown,  R.C.  and  Kolb,  C.E., 

"Kinetics  of  High  Temperature  B/0/H/C  Chemistry,"  Combustion  and  Flame 
M,  43  (1991). 


3-26 


3-4  Brown,  R.C.,  Kolb,  C.E.,  Yetter,  R.A.,  Dryer,  F.L.  and  Rabitz,  H.R., 
"Development  of  a  Boron  Assisted  Combustion  Model  with  Sensitivity 
Analysis,"  Aerodyne  Research,  Inc.  Final  Report  No.  ARI-RR-580  (1987). 

3-5  Page,  M. ,  "Multireference  Configuration  Interaction  Study  of  the 
Reaction  H2  +  BO  -*■  HBO,"  J.  Phys.  Chem.  93,  3693  (1989). 

3-6  JANAF  Thermochemical  Tables,  Third  Edition,  Parts  I  and  II,  J.  of 
Physical  and  Chemical  Reference  Data,  Volume  14  (1985). 

3-7  Wedler,  G.,  Chemisorption:  An  Experimental  Approach.  Butterworths , 
(1976). 

3-8  Turns,  S.R. ,  Holl,  J.T. ,  Solomon,  A.S.P.  and  Faeth,  G.M.,  "Gasification 
of  Boron  Oxide  Drops  in  Combustion  Gases,"  Combust.  Sci.  and  Tech.  43. 
287  (1985). 

3-9  MacKenzie,  J.D.,  "Structure  of  Liquid  Boron  Trioxide, "  J.  Chem.  Phys. 
22,  605  (1958). 

3-10  Silver,  A.H.  and  Bray,  P.J.,  "Nuclear  Magnetic  Resonance  Absorption  in 
Glass  I.  Nuclear  Quandruple  Effects  in  Boron  Ox'de,  Soda  Boric  Oxide 
and  Borosilicate  Glasses,"  J.  Chem.  Phys.  2i»  984  (1958). 

3-11  Silver,  A.H.,  "Nuclear  Magnetic  Resonance  in  B2O3-H2O  Glasses  and  Boric 
Acids,"  J.  Chem.  Phys.  22,  959  (1960). 

3-12  Krogh-Moe,  J.  ,  "interpretation  of  the  Infrared  Spectra  of  Boron  Oxide 
and  Alkali  Borate  Glasses,"  Phys.  and  Chem.  Glasses  £,  46  (1965). 

3-13  Mozzi,  R.L.  and  Warren,  B.E.,  "The  Structure  of  Boron  Oxide,"  J.  Appl. 
Cryst.  3,  251  (1970). 

3-14  Jellison,  Jr.,  G.E.,  Panek,  L.W.,  Bray,  P.J.  and  Rouse,  Jr.,  G.B., 
"Determination  of  Structure  and  Bonding  in  Vireous  B2O3  by  Means  of 
B10,  B11,  and  017  NMR,"  J.  Chem.  Phys.  ££,  802  (1988). 

3-15  Miyake,  M. ,  Suzuki,  T.,  Morikawa,  H.,  Takayi,  and  Marumo,  F., 

"Structural  Analysis  of  Molten  B2O3,"  J.  Chem.  Soc.  Faraday  Trans.  1, 
M,  1925  (1984). 

3-16  Elliott,  S.R.,  "A  Continuous  Random  Network  Approach  to  the  Structure 
of  Vitreous  Boron  Trioxide,"  Phil.  Mag.  B  22,  435  (1978). 

3-17  Soules,  T.F. ,  "A  Molecular  Dynamic  Calculation  of  the  Structure  of  B2O3 
Glass,"  J.  Chem.  Phys.  72,  4032  (1980). 


3-27 


3-18  Glassman,  I.,  Williams,  F.A.  and  Antaki,  P.,  "A  Physical  and  Chemical 
Interpretation  of  Boron  Particle  Combustion,"  Twentieth  Symposium 
(International')  on  Combustion,  the  Combustion  Institute,  Pittsburgh,  PA 
(1984). 

3-19  Hanley,  L.,  Whitten,  J.L.  and  Anderson,  S.L. ,  "Collision- Induced 

Dissociation  and  ab  Initio  Studies  of  Boron  Cluster  Ions:  Determination 
of  Structures  and  Stabilities,"  J.  Phys.  Chem.  92,  5803  (1988). 

3-20  Hanley,  L.  and  Anderson,  S.L.,  "Oxidation  of  Small  Boron  Clusters  Ions 
(B+]_i3)  by  Oxygen,"  J.  Chem.  Phys.  89,  2848  (1988). 

3-21  Ruatta,  S.A.,  Hanley,  L.  and  Anderson,  S.L. ,  "Dynamics  of  Boron  Cluster 
Ion  Reactions  with  Deuterium:  Adduct  Formation  and  Decay,"  J.  Chem. 
Phys.  91,  226  (1989). 

3-22  Hintz,  P.A.,  Ruatta,  S.A.  and  Anderson,  S.L.,  "interaction  of  Boron 
Cluster  Ions  with  Water:  Single  Collision  Dynamics  and  Sequential 
Etching,"  J.  Chem.  Phys.  92,  292  (1990). 

3-23  Hintz,  P.A.  Ruatta,  S.A.  and  Anderson,  S.L.,  "Cluster  Based  Reaction 

Probabilities  for  Boron  with  Oxygen,  Hydrogen,  Water,  Nitrogen,  Nitrous 
Oxide,  Carbon  Dioxide,  Carbon  Monoxide,  Methane,  Teltraf luoromethane 
and  Silane,"  Interim  Report  SBIL-A89,  State  University  of  New  York, 
Chemistry  Department  (1989). 

3-24  Gomez,  A.,  Zvuloni,  R.  and  Rosner,  D.E.,  "Flow  Reactor  Studies  of  Boron 
Vaporization  and  Kinetically-Controlled  Oxidation  Using  Microwave- 
Induced  Plasma  Emission  Spectroscopy:  Preliminary  Results,"  presented 
at  the  1986  Technical  Meeting.  The  Combustion  Institute  Eastern  States 
Section,  Dec.  12,  1986,  San  Juan,  Peurto  Rico,  and  at  the  Joint  Meeting 
of  the  French  and  Italian  Section  of  the  Combustion  Institute,  June  16, 
1987,  Amalfi,  Italy. 

3-25  Zvuloni,  R. ,  Gomez,  A.  and  Rosner,  D.E.,  "Direct  Measurements  of  the 
High  Temperature  Kinetics  of  Solid  Boron  Gasification  by  Its  Higher 
Oxide  8203(g):  Chemical  Propulsion  Implications,”  J.  Propulsion  and 
Power,  in  press  (1989). 

3-26  Zvuloni,  R.  and  Rosner,  D.E.,  "High  Temperature  Gasification  Rate  of 

Solid  Boron  In  Mixtures  of  its  Highest  Oxide  8203(g),  and  Water  Vapor," 
AIAA  J.  Propulsion  and  Power,  Submitted  1989. 


3-28 


4.0  MODEL  RESULTS 


4 . 1  B2Q3X&}  Gasification 

The  spherical  particle  combustion  model  described  in  Section  2.0  was  used 
to  simulate  the  gasification  of  a  liquid  boron  oxide  droplet  injected  into 
high  temperature,  hydrocarbon  combustion  environment. 

In  the  flat-flame  burner  experiments  of  a  single  suspended  droplet,^-1 
droplet  diameters  were  on  the  order  of  1000  pm.  In  contrast,  much  smaller 
(20  -  100  pm)  boron  particles  coated  with  liquid  boron  oxide  have  been  used  to 
observe  ignition  characteristics .^“2  jn  the  present  model  simulations,  there¬ 
fore,  the  droplet  diameter  was  varied  between  50  -  1000  pm  to  examine  the 
effects  of  the  droplet  size.  In  addition,  ambient  temperatures  between  1600  K 
and  2200  K  were  considered.  This  temperature  range  covers  the  observed 
ignition  temperature  of  particulate  boron.  The  ambient  gas  mixture 
composition  at  three  temperatures  is  presented  in  Table  13.  This  data  was 
obtained  from  equilibrium  calculations  for  a  JP4/air  mixture  with  an 
equivalence  ratio  of  0.5. 


Table  13  -  The  Ambient  Mixture  Compositions  in  Mole  Fractions  Obtained  from 
Equilibrium  Calculation  of  a  JP-4/Air  Mixture  with  an  Equivalence 
Ratio  of  0.5 


Temperature 

1600  K 

1800  K 

2000  K 

0 

0.0005% 

0.004% 

0.02% 

OH 

0.019% 

0.05% 

0.14% 

°2 

10% 

10% 

10% 

h2o 

6.8% 

6.8% 

6.7% 

C02 

6.8% 

6.8% 

6.8% 

Temporal  variations  of  the  gasification  rate  and  surface  species 
concentrations  for  a  50  pm  boron  oxide  droplet  and  an  ambient  temperature  of 
1800  K  are  presented  in  Figure  4.1.  Figure  4.1  shows  that  the  initial 
transient  period  is  relatively  short;  less  than  0.05%  of  the  particle's 
lifetime  tf.  Additional  model  simulation  results  confirmed  this  observation 
for  the  range  of  conditions  treated  here.  Thus,  only  the  quasi-steady  state 
behavior  will  be  discussed  hereafter  and  all  model  results  for  gasification 
rates,  surface  chemical  species  concentrations,  gas  phase  species 
concentrations,  temperatures,  etc.  are  the  quasi-steady  state  values  unless 
otherwise  noted. 


t/tf 

Figure  4.1.  Temporal  Variation  of  the  6203(2)  Gasification  Rate  and  Percent 
Surface  Coverage  By  Adsorption  Products  for  a  50  pm  Particle 
at  an  Ambient  Temperature  of  1800  K.  tf  denotes  the  droplet 
lifetime. 
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Figures  4.2  and  4.3  show  the  gasification  rates  due  to  vaporation  and 
surface  reactions  with  H2O/OH  versus  droplet  diameter  for  four  ambient 
temperatures.  These  plots  can  be  used  to  identify  rate  controlling  steps 
since  the  gasification  rate  stays  constant  in  the  diffusion  controlled  regime 
and  is  linearly  proportional  to  the  particle  diameter  in  the  kinetically 
controlled  regime.  At  ambient  temperatures  of  1600  K  and  1800  K,  Figure  4.2 
shows  that  gasification  by  surface  chemical  reactions  is  kinetically 
controlled  for  particles  smaller  than  200  pm.  As  the  particle's  diameter 
increases,  the  diffusion  process  slows  down  and  eventually  becomes  the  rate 
controlling  step.  In  contrast,  gasification  by  the  physical  evaporation  is 
always  diffusion  controlled. 

The  gasification  rates  as  a  function  of  particle  diameter  for  ambient 
temperatures  of  T  =  2000  K  and  T  =  2200  K  are  shown  in  Figure  4.3.  In 
contrast  to  the  low  ambient  temperature  cases  (T  =  1600  and  1800  K)  presented 
in  Figure  4.2,  at  these  higher  temperatures  the  evaporation  gasification  rate 
increases  sharply  for  particle  diameters  less  than  200  pm.  This  is  due  to  the 
strong  effect  of  the  surface  temperature  on  evaporation  and  not  to  changes  in 
the  rate-controlling  steps.  Analogous  to  the  low  temperature  cases, 
gasification  by  surface  reactions  is  kinetically  limited 
while  the  gasification  by  the  evaporation  is  diffusion  limited  for  small 
particles.  Therefore,  as  the  particle  diameter  decreases,  the  heat  loss  per 
unit  area  due  to  the  surface  reactions  stays  the  same  but  the  diffusional  heat 
transfer  from  the  ambient  to  the  particle  increases.  As  a  result,  the  surface 
temperature  increases  as  the  particle  diameter  decreases.  This  surface 
temperature  increase  is  large  enough  to  increase  the  physical  evaporation  rate 
substantially  under  high  ambient  temperature  conditions. 

Figures  4.2  and  4.3  also  show  that  as  the  ambient  temperature  increases, 
gasification  by  the  physical  evaporation  is  accelerated  faster  than 
gasification  by  the  surface  reactions.  This  is  due  to  differences  in  the 
associated  activation  energies. 
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Figure  4.2.  Calculated  8203(H)  Gasification  Rates  Due  to  Surface  Reaction 
and  Physical  Evaporation  as  a  Function  of  Particle  Diameter 
at  Ambient  Temperatures  of  (a)  1600  K  and  (b)  1800  K 
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(b) 


Figure  4. 


.  Calculated  6203(H)  Gasification  Rates  Due  to  Surface  Reaction 
and  Physical  Evaporation  as  a  Function  of  Particle  Diameter 
at  Ambient  Temperatures  of  (a)  2000  K  and  (b)  2200  K 
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Figures  4.2  and  4.3  also  demonstrate  that  gasification  by  the  surface 
reactions  dominates  over  over  the  physical  evaporation.  This  is  especially 
true  for  larger  particles  where  both  mechanisms  are  diffusion  controlled. 
However,  gasification  by  the  surface  reactions  becomes  slower  as  the  particle 
diameter  decreases.  Figure  4.3  shows  that  physical  evaporation  is  more 
important  than  surface  reactions  in  gasifying  the  small  droplets  (<  100  pm)  at 
temperatures  above  2000  K.  Consequently,  experimental  results  obtained  for 
1000  uni  particles  by  Turns  et  al.^"*  should  be  interpreted  very  carefully, 
because  their  conclusions  about  the  relative  importance  of  surface  reactions 
may  not  be  directly  applicable  to  smaller  particles. 

The  overall  gasification  rate  are  obtained  by  simply  adding  the 
individual  gasification  rates  due  to  surface  reactions  and  physical 
evaporation.  Figure  4.4  shows  the  overall  gasification  rate  as  a  function  of 
particle  diameter  for  ambient  temperatures  of  1600  K,  1800  K,  2000  K  and 
2200  K.  For  large  particles  (>  200  pm) ,  gasification  is  diffusion  limited  and 
the  overall  gasification  rates  are  independent  of  particle  diameter.  For 
small  particles  (<  200  pm),  gasification  is  kinetically  limited  for  the  low 
ambient  temperatures  (T  =  1600  K  and  1800  K) ,  but  deviates  considerably  from 
the  dl  behavior  for  the  high  ambient  temperatures  (T  =  2000  K  and  2200  K) . 

This  deviation  from  d*  behavior  is  caused  by  the  increased  contribution  from 
physical  evaporation  as  shown  in  Figure  4.3. 

Gas  phase  chemical  species  and  temperature  profiles  are  shown  in 
Figure  4.5  for  particle  diameters  of  50  pm  and  500  pm  and  an  ambient 
temperature  of  1800  K.  The  boron  containing  species,  B0£(g),  8203(g),  and 
HOBO(g)  are  mainly  produced  from  liquid  boron  oxide  surface  reactions  and  thus 
have  their  maximum  values  near  the  surface.  B02(g)  shows  the  largest  decrease 
with  distance  from  the  surface,  and  the  fall  off  is  particularly  sharp  for  the 
500  pm  diameter  particle  (Figure  4-5b).  This  is  because  for  large  particles, 
diffusion  mass  transfer  is  slow  relative  to  the  homogeneous  gas  phase 
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Figure  4.4.  Calculated  Overall  8203(11)  Gasification  Rates  as  a  Function  of 
Particle  Diameter  at  Ambient  Temperatures  of  1600  K,  1800  K, 

2000  K,  and  2200  K. 

chemistry  that  converts  B02(g)  to  HOBO  and  8203(g).  The  higher  concentrations 
of  B02(g)  and  HOBO(g)  at  the  surface  of  the  500  ym  diameter  particle  are  also 
due  to  slower  diffusion  mass  transfer  rates.  In  contrast,  the  8203(g) 
concentration  at  the  surface  of  the  500  ym  diameter  particle  is  found  to  be 
lower  than  that  of  the  50  ym  diameter  particle  because  of  the  lower  surface 
temperature. 

As  listed  in  Section  3.0,  the  surface  reaction  mechanisms  used  in  the 
present  work  consist  of  adsorption,  reaction  between  adsorbed  species  and 
desorption.  The  reaction  fluxes  from  each  elementary  step  are  presented  in 
Figures  4.6,  4.7,  and  4.8.  These  reaction  fluxes  can  be  used  to  identify 
major  pathways. 
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Figure  A. 5.  Radical  Species  Mass  Fraction  and  Temperature  Profiles  for  an 

Ambient  Temperature  of  1800  K  and  Droplet  Diameter  of  (a)  50  pm 
and  (b)  500  pm. 
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Adsorption  and  desorption  fluxes  as  a  function  of  particle  diameter  are 
shown  in  Figures  4.6  and  4.7.  As  the  droplet  diameter  increases,  the 
adsorption  fluxes  of  HOBO(g)  and  the  desorption  fluxes  of  >B-OH(c)  increase 
rapidly  in  the  kinetically  controlled  regime  and  increase  only  slightly  in  the 
diffusion  controlled  regime.  They  also  increase  rapidly  as  the  ambient 
temperature  increases.  Even  though  the  adsorption  of  HOBO(g)  has  the  largest 
adsorption  flux,  the  desorption  of  >B-0H(s)  is  also  fast  enough  to  return  HOBO 
molecules  back  to  the  gas  phase.  The  adsorption  of  H20(g)  appears  to  have  the 
largest  flux  next  to  that  of  HOBO.  As  shown  in  Table  4,  the  adsorption  of  one 
H20(g)  molecule  creates  two  surface  species  >B-OH(c)  which  desorb  to  produce 
two  molecules  of  HOBO(g).  The  adsorption  flux  of  0(g)  atom  and  B02(g),  as 
well  as  the  desorption  flux  of  >B0(c)  is  small  for  the  given  ambient 
conditions  in  which  the  concentration  0(g)  is  much  smaller  than  that  of  0H(g) 
or  H20(g).  Because  of  the  abundance  of  0H(g)  and  H20(g)  in  the  ambient 
environment,  the  desorption  of  >B-0H  is  the  dominant  desorption  pathway  in 
gasifying  liquid  boron  oxides.  This  is  shown  in  Figure  4.7. 

The  reaction  fluxes  for  second  order  reactions  between  two  surface 
species  are  shown  in  Figure  4.8.  They  were  found  to  be  one  order  of  magnitude 
smaller  than  the  reaction  fluxes  for  adsorption  or  desorption.  Furthermore, 
the  calculated  surface  coverage  by  the  adsorbed  species  are  very  small  (<  IX) 
as  shown  in  Figure  4.9.  This  is  because  desorption  is  sufficient  (relative  to 
adsorption)  to  keep  the  surface  coverage  low. 

Sensitivity  analysis  is  a  useful  tool  in  identifying  important  input 
parameters  in  determining  physically  important  output  parameters.  In  the 
present  work,  a  parameter  sensitivity  coefficient  of  the  gasification  rate 
with  respect  to  surface  kinetic  rate  constants  was  used  to  measure  the 
response  of  the  gasification  rate  to  variation  in  the  surface  reaction  rate 
parameters.  This  analysis  showed  that  the  rate  parameters  associated  with 
H20(g)  and  0H(g)  adsorption,  >B-0H(s)  desorption  and  the  second  order  reaction 
>B-0H  +  >B-0H  are  important  in  determining  gasification  rate.  Figure  4.10 
displays  sensitivity  coefficient  profiles  for  these  important  rate  parameters 
as  a  function  of  droplet  diameter  for  an  ambient  temperature  1800  K. 
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Figure  4.6.  Adsorption  Fluxes  as  a  Function  of  Droplet  Diameters  for  Ambient 
Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K  and 
(d)  2200  K  (Continued) 
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Figure  4.7.  Desorption  Fluxes  as  a  Function  of  Droplet  Diameters  for  Ambient 
Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K  and 
(d)  2200  K. 
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Figure  4.7.  Desorption  Fluxes  as  a  Function  of  Droplet  Diameters  for  Ambient 
Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K  and 
(d)  2200  K  (Continued) 
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Figure  4.8.  Surface  Reaction  Flux  as  a  Function  of  Droplet  Diameter  for 

Ambient  Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K  and 
(d)  2200  K 
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Figure  4.8. 
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Surface  Reaction  Flux  as  a  Function  of  Droplet  Diameter  for 
Ambient  Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K  and 
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Figure  4.9.  Percentage  Surface  Coverage  as  a  Function  of  Droplet  Diameter 
for  Ambient  Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K 
and  (d)  2200  K 
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Figure  4.9.  Percentage  Surface  Coverage  as  a  Function  of  Droplet  Diameter 
for  Ambient  Temperatures  of  (a)  1600  K,  (b)  1800  K,  (c)  2000  K 
and  (d)  2200  K  (Continued) 
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Figure  A. 10.  Sensitivity  Coefficient  Profiles  as  a  Function  of  Droplet 
Diameter  for  an  Ambient  Temperature  of  1800  K 

4.2  Hjgb  Temper^turq  p^ticulate  BQtpn 

The  section  presents  model  results  for  the  oxidation  of  200  ym  solid 
boron  particle  in  a  hydrocarbon  combustion  environment  using  the  surface 
reactions  and  rate  parameters  described  in  Section  3.5.  The  same  ambient 
mixture  composition  (Table  13  in  Section  4.1)  used  in  the  model  simulations  of 
liquid  B2O3  gasification  was  also  used  to  treat  particulate  boron  oxidation. 
Model  results  are  presented  here  for  two  cases  corresponding  to  ambient 
temperatures  of  1400  K  and  2000  K. 

For  the  simulations  described  here,  the  particle  radius  was  kept  fixed  in 
order  to  explicitly  examine  the  steady-state  behavior.  The  burning  rates  for 
the  two  cases  studied  are  presented  in  Figure  4-11.  The  transient  period 
prior  to  steady-state  conditions  is  clearly  evident.  The  final  time,  tf,  in 
this  figure  was  chosen  as  10  sec.  It  is  seen  that  for  both  Tamb  =  1400  K  and 
2000  K,  a  steady-state  is  reached  within  12  msec.  Figure  4-11  also  indicates 


4-18 


that  the  burning  rate  is  only  weakly  dependent  on  the  ambient  temperature.  In 
particular,  a  600  K  increase  in  the  ambient  temperature  only  increases  the 
burning  rate  constant  by  ca.  11%. 

Temperature  profiles  are  shown  in  Figure  4-12.  These  profiles  indicate  a 
temperature  at  the  particle  surface  that  is  between  600  K  and  700  K  higher 
than  the  ambient.  This  difference  is  primarily  dictated  by  surface  reaction 
enthalpies  and,  consequently,  subject  to  the  current  uncertainties  in 
thermochemical  parameters.  The  temperature  profiles  for  both  cases  treated 
decrease  monotonically  with  radial  distance  from  surface  and  fall  to  within 
100  K  of  the  ambient  temperature  within  nine  particle  radii. 

Radial  profiles  of  the  species  mass  fractions  for  Tamb  =  1400  K  are  shown 
in  Figure  4-13.  It  is  readily  seen  that  the  mass  fractions  for  gas  phase 
oxidizers  decrease  near  the  particle  surface.  Although  gas  phase  chemistry 
may  play  some  role  in  this  behavior,  the  decrease  is  primarily  attributable  to 
reactions  at  the  surface.  Consequently,  the  extent  of  this  decrease  can  be 


Figure  4.11.  Calculated  B(s)  Burning  Rate  Constant  for  an  Ambient 
Temperature  of  1400  K  and  2000  K 
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used  as  a  qualitative  criterion  for  judging  the  relative  oxidizing  strength  of 
the  various  gas  phase  reactants.  On  this  basis,  the  model  results  agree 
qualitatively  with  observed  experimental  reaction  probabilities.  In 
particular,  0(g),  0H(g)  and  8203(g)  are  seen  to  be  relatively  strong 
particulate  boron  oxidizing  agents,  while  002(g),  02(g)  and  H20(g)  are 
relatively  moderate.  H0B0(g)  was  clearly  found  to  be  the  least  effective  gas 
phase  reactant,  however,  there  is  little  experimental  rate  data  for 
B(s)/H0B0(g)  reactions  currently  available  and  the  estimated  rate  parameters 
given  in  Section  3.0  are  highly  uncertain.  In  addition,  because  the  H0B0(g) 
mass  fraction  is  significantly  larger  then  the  mass  fractions  of  other 
potential  reactants,  it  may  be  an  imrnrtant  oxidizing  species  even  if  its  rate 
constant  is  relatively  small. 

The  mass  fractions  for  surface  reaction  products  B02(g),  00(g),  B0(g), 
8203(g),  HBO(g),  H(g)  and  H2(g)  all  show  pronounced  peaks  at  the  surface.  The 
species  containing  boron,  however,  decrease  rather  sharply  and,  with  the 
exception  of  B02(g),  are  negligible  beyond  one  particle  radius  (200  pm)  from 
the  surface.  This  is  the  consequence  of  homogenous  gas  phase  chemistry  which 
converts  boron  sub-oxides  and  oxyhydrides  to  the  primary  boron  oxidation 
products  H0B0(g)  and  8203(g).  Here  there  is  good  agreement  with  previous 
homogeneous  oxidation  model  results  which  indicated  that  surface  products  were 
quickly  converted  to  H0B0(g)  and  8203(g).  Similarly,  the  indication  in 
Figure  4-13  that  H0B0(g)  tends  to  peak  closer  to  the  surface  than  8203(g)  is  a 
consequence  of  the  gas  phase  chemistry  model  which  was  shown  to  predict  that 
H0B0(g)  was  formed  more  rapidly  by  gas  phase  reactions  than  8203(g). 

Radial  profiles  of  the  species  mass  fractions  for  Tamb  =  2000  K  are  shown 
in  Figure  4-14.  In  general,  the  mass  fraction  profiles  exhibit  the  same 
trends  as  described  above  for  the  case  corresponding  to  Tamb  *  1400  K.  The 
most  notable  difference  is  the  decrease  in  the  H0B0(g)  mass  fraction  and  a 
corresponding  increase  in  the  8203(g)  mass  fraction.  This  behavior  is  a 
direct  consequence  of  the  gas  phase  chemistry  which  was  found  in  earlier 
modeling  work  to  show  a  shift  towards  8203(g)  with  higher  ambient 
temperatures. 
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Figure  4.13.  Species  Mass  Fraction  Radial  Profiles  for  the  Oxidation 
of  a  200  ym  Sperical  Boron  Particle  in  a  JP4/air  Mixture 
at  an  Ambient  Temperature  of  1400  K 
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Figure  4.13.  Species  Mass  Fraction  Radial  Profiles  for  the  Oxidation 
of  a  200  ym  Sperical  Boron  Particle  in  a  JP4/air  Mixture 
at  an  Ambient  Temperature  of  1400  K  (Continued) 
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Figure  4.13.  Species  Mass  Fraction  Radial  Profiles  for  the  Oxidation 
of  a  200  pm  Sperical  Boron  Particle  in  a  JP4/air  Mixture 
at  an  Ambient  Temperature  of  1400  K  (Continued) 
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Figure  4.14.  Species  Mass  Fraction  Radial  Profiles  for  the  Oxidation  of  a 
200  ym  Spherical  Boron  Particle  in  a  JP4/air  Mixture 
at  an  ambient  temperature  of  2000  K 


4-25 


mass  fraction 


.005 


Figure  4.14.  Species  Mass  Fraction  Radial  Profiles  for  the  Oxidation  of  a 
200  pm  Spherical  Boron  Particle  in  a  JP4/air  Mixture 
at  an  ambient  temperature  of  2000  K 
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Figure  A. 14.  Species  Mass  Fraction  Radial  Profiles  for  the  Oxidation  of  a 
200  pm  Spherical  Boron  Particle  in  a  JP4/air  Mixture 
at  an  ambient  temperature  of  2000  K  (Continued) 

Overall,  therefore,  the  model  results  indicate  that  the  species 
distribution  near  the  particle  surface  is  predictably  governed  by 
heterogeneous  surface  reactions.  However,  homogeneous  gas  phase  chemistry 
rapidly  converts  the  species  formed  through  surface  reactions  to  the  primary 
boron  oxidation  products  of  H0B0(g)  and  8203(g).  At  radial  distances  greater 
than  ca.  one  particle  radius  from  the  surface,  the  species  distribution  and 
temperature  is  governed  by  gas  phase  oxidation  reactions.  Because  the 
products  of  surface  reactions  undergo  gas  phase  reactions  so  quickly,  the  gas 
phase  chemistry  is  largely  independent  of  the  surface  chemistry.  The  feature 
is  consistent  with  earlier  modeling  results  for  the  gas  phase  oxidation 
model^-^  and  a  continuously  stirred  reactor  model^"^  for  coupling  between 
surface  reactions  and  gas  phase  chemistry. 
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